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CALCULATION OF SUPERSONIC STREAM PARAMETERS OF A REAL GAS 
FROM MEASURABLE QUANTITIES USING FORTRAN IV ROUTINES 

by Robert C. Johnson 
Lewis Research Center 

SUMMARY 

Two sets of FORTRAN IV routines are presented that calculate flow and thermody- 
namic properties of a supersonic stream from measurable quantities. These flow and 
thermodynamic properties include velocity, density, enthalpy, entropy, and isentropic 
exponent, both in the free stream and behind the normal shock. Two sets of measure- 
ments are considered. For one set, these measurements are the stagnation pressure, 
the stagnation temperature, and the pressure on the surface of a static -pressure wedge. 
For the other set, they are the pressure and temperature in a plenum upstream of the 
supersonic nozzle and the stagnation pressure at the exit of this nozzle. These routines 
apply to any undissociated gas whose real-gas properties are known. As examples, the 
routines are specifically applied to air, nitrogen, oxygen, normal hydrogen, parahydro- 
gen, helium, argon, steam, methane, and natural gas. 


INTRODUCTION 

When the properties of a supersonic stream are to be determined from measured 
quantities, there are a number of sets of quantities that can be measured. For example, 
one set, which is useful when local stream properties are desired, could be the local 
stagnation pressure and temperature and the static pressure on the surface of a cone or 
wedge. (The term stagnation refers to total conditions downstream of the normal shock. ) 
If the supersonic flow issues from a nozzle connected to an upstream plenum, another set 
of measurements is useful when the flow from the plenum to the nozzle exit can be con- 
sidered isentropic and one dimensional (e. g. , shock -free flow outside of the boundary 
layer). These measurements could be the pressure and temperature within the plenum 
and the stagnation pressure at the nozzle exit. From these measurements, such stream 
properties as Mach number, velocity, density, and temperature can be determined. How- 
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ever, since these stream properties are implicity rather than explicitly related to the 
measured quantities, iterative procedures must be used to make these determinations. 
For the case of a perfect gas, tables and graphs exist that simplify these calculations 
(e.g. , ref. 1). In this report, a perfect gas is defined as one whose compressibility 
factor equals unity and whose specific heat is invariant. This perfect gas is to be distin- 
guished from an ideal gas, which like a perfect gas is defined to have a compressibility 
factor of unity but, unlike a perfect gas, is defined to have a temperature -dependent 
specific heat. In the absence of dissociation, all gases approach this ideal-gas condition 
as the pressure is reduced. A nonperfect gas is called a real gas. For a gas such as 
air at room temperature and at pressures less than a few atmospheres, the perfect-gas 
assumption is adequate, at least for moderate Mach numbers. For gases at high pres- 
sures, and/or low temperatures, the perfect-gas assumption is not sufficiently valid for 
accurate computations, and real-gas computations should be made. Such calculations 
were made in reference 2, where high-pressure nitrogen was used to generate a Mach 3 
stream . 

In this report, the FORTRAN IV computer routines that were used to make the cal- 
culations in reference 2 are presented. In one set of routines, the independent variables 
are (1) the pressure and temperature such as would be measured by a stagnation pressure 
probe and a stagnation temperature probe and (2) the static pressure that would be meas- 
ured on the surface of a wedge. In the other set of routines, the independent variables 
are (1) the total pressure and total temperature that would be measured in the subsonic 
plenum upstream of the supersonic nozzle and (2) the pressure that would be measured 
by a stagnation pressure probe at the exit of the nozzle. These routines apply to the 
following gases: air, nitrogen, oxygen, normal hydrogen, parahydrogen, helium, argon, 
steam, methane, and natural gas. The sections of the routines that have to do with the 
nature of the gas can also be used in conjunction with the routines given in reference 3. 
The combined set of routines can then be used to calculate the mass flow rate of these 
real gases through subsonic or sonic flow nozzles. 

The pressure and temperature ranges for which these routines give accurate results 
depend on the pressure and temperature ranges of the state equations used to describe 
the behavior of the various gases. These ranges are given in table I for each gas. Gen- 
erally, these routines are accurate at very low pressures, as long as the flow is in the 
continuum regime and no gas dissociation is present. Continuum flow exists when the 
mean free path is small compared to the probe dimensions. Dissociation is most liable 
to occur at high temperatures and low pressures. The lower limit of pressure is gener- 
ally set by the criterion of continuum flow; the other limits of pressure and temperature 
are set by the range of validity of the state equations. In a supersonic stream, there can 
be a large pressure and temperature variation of the gas as it accelerates from a plenum, 
where the gas is at rest, to the exit of the supersonic nozzle. These pressures and tem- 
peratures should all be within the pressure and temperature ranges prescribed in table I. 
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TABLE I- - PRESSURE AND TEMPERATURE RANGES 


FOR THE VARIOUS GASES 


Gas 

Pressure range, 
pascal 

Temperature range, 
K 

Air 

100 to lOOxlO 5 

150 to 1500 

Nitrogen, range I 

100 to 300X10 5 

60 to 400 

Nitrogen, range II 

100 to 100x10 s 

170 to 3000 

Oxygen, range I 

100 to 300xl0 5 

60 to 400 

Oxygen, range II 

100 to lOOxlO 5 

180 to 3000 

Normal hydrogen 

100 to 100x10 s 

70 to 600 

Parahydrogen 

flOO to 300X10 5 

13 to 100 


llOO to 100x10 s 

13 to 400 

Helium 

100 to 300X10 5 

5.4 to 400 

Argon 

100 to 1000x10 s 

30 to 400 

Steam 

100 to 1000X10 5 

273 to 1600 

Methane 

100 to 400x10 s 

70 to 600 

Natural gas 

100 to lOOxlO 5 

200 to 400 


In addition to such quantities as the free-stream velocity and density, additional 
thermodynamic quantities such as enthalpy, entropy, specific heat at constant pressure, 
specific -heat ratio, and isentropic exponent are calculated not only for free-stream con- 
ditions but also for plenum conditions, stagnation conditions behind the normal shock, 
free-stream conditions behind the normal shock, and (when applicable) free-stream con- 
ditions behind the oblique shock. 

Both sets of routines are in two sections: one section is concerned with the iteration 
procedures, and the other section involves the specific behavior of the gas. The itera- 
tion routines that apply when the independent variables are the stagnation pressure, the 
stagnation temperature, and the static pressure at the surface of a wedge are described 
and presented in appendix B. The iteration routines that apply when the independent 
variables are the plenum pressure, the plenum temperature, and the stagnation pressure 
are described and presented in appendix C. For each of the gases, appendix D describes 
and presents the routines that have to do with the specific behavior of the gas. In addi- 
tion, appendix D discusses the state equations that are involved in these routines. The 
routines in appendix D are used in conjunction with the routines in either appendix B or 
C and can also be used in conjunction with the iteration routines in reference 3 to calcu- 
late the mass flow rate of these gases through subsonic or sonic flow nozzles. Appen- 
dix E includes sample calculations for a supersonic air stream. 

All symbols are defined in appendix A. The International System of Units (SI) is 
used throughout this report, except that appendix E illustrates the conversion from U. S. 
customary units to SI units in a practical example. 
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ANALYSIS 


Basic Relations 


The compressibility factor for all the gases involved in this report is given as a 
function of density and temperature. Therefore, this analysis requires that such thermo- 
dynamic functions as entropy, enthalpy, and specific heat also be expressed as a function 
of density and temperature. This involves the following compressibility -factor functions: 


Zj (p, T) = Z = — — 

1 pRT 

Z n (p,T)= Z + t(^) =~(—) 

11 \dT/ p Rp \0T l p 

Z m (p.T,= Z + pg-) T -£(*)_ 
Z IV (p,T)= f P { Z n - 1)& 

Jo p 

Z V <P,T)= f P { Zjj-Zj)^ 

*r\ P 


^ Vi^P? T) - 


. T / az iv\ 

v 3t A 


(1) 

( 2 ) 

(3) 

(4) 

(5) 

( 6 ) 


For an ideal gas, Z V Z H’ and Zj n equal 1 and Zjy, Zy, and Zyj equal zero. In 
addition, for any real gas, Zj, Zjj, and Z^j have to be positive . 

In addition, the following functions of the ideal-gas specific heat C v are involved: 


^(T) = 



(7) 


g w (T) = f C *’ id g. + ln-P- ( 8 ) 

11 •/ R T R RT 

^ m (T) = f -Jf> id dT = ^ - T (9) 

J R R 
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In terms of the functions defined by equations (1) to (9), the following thermodynamic 
functions can now be expressed in terms of density and temperature: 


- = ZjpT 

R 1 

(10) 

|= 4 n -lnp-Z IV 

(ii) 

f=«m +T < z i- z v> 

(12) 

¥ 

(13) 

5 > = ^ + Z 1 

R R Zjjj 

(14) 

n 

Ol o 

< TO 

(15) 



where k is the isentropic exponent, and 

2 

2 — = kZ T T (17) 

R 1 

Except for some symbol changes, equations (1) to (17) can be found in reference 4. The 
compressibility factor and the ideal-gas specific heat for the gases involved in this report 
are discussed in appendix D. 

Now that the thermodynamic functions given by equations (10) to (17) are in terms of 
density and temperature, the procedures for calculating supersonic stream properties 
can be discussed. As indicated in the INTRODUCTION, there are two cases considered, 
and each one requires a different procedure. In the first case, measurements are made 
by a stagnation pressure probe, a stagnation temperature probe, and a wedge-type static- 
pressure probe. The quantities represented by these measurements are the total pres- 
sure and total temperature downstream of the normal shock (assuming that the recovery 
correction is applied to the temperature measurement) and the static pressure downstream 
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of the oblique shock attached to the wedge. In the second case, the measurements are 
the pressure and temperature in the plenum upstream of a supersonic nozzle and the 
pressure indicated by a stagnation pressure probe at the nozzle exit. The quantities rep- 
resented by these measurements are the free-stream total pressure and total temperature 
and the total pressure downstream of the normal shock. 

These cases will be discussed separately. A word, first, about the notation. The 
subscript 1 refers to free-stream conditions, the subscript 2 refers to conditions down- 
stream of the normal shock, and the subscript 3 refers to conditions downstream of 
the oblique shock that is attached to the wedge. The second subscript, t, when appended, 
refers to total conditions (i. e. , conditions that would exist if the gas were brought to rest 
isentropically) . Thus, a symbol like Hg refers to the enthalpy downstream of the nor- 
mal shock, and a symbol like Hg ^ refers to the total enthalpy downstream of the nor- 
mal shock. The computational procedures for both of these cases involve solving nonlin- 
ear equations or sets of nonlinear equations. The Newton-Raphson iteration procedure is 
used in all cases. 

Case I - given Pg t , Tg Pg, and 5. - Since the thermodynamic equations involve 

density and temperature rather than pressure and temperature, the following equation is 
implicity solved for pg t : 

p 2,t = p * p 2,t’ T 2,P ( 18 > 

Next, there are two equations that relate total conditions to static conditions downstream 
of the normal shock. These are 


s 0>2,t> T 2,t>- S( f2> T 2> 

(19) 

’2,t’ T 2,t* “ H ^2, T 2* * " U 2 

(20) 


Then, the energy, momentum, and continuity equations are written for the flow across 
the normal shock. These are 


H(p r Tj) + i Uj = H (p 2 , T 2 ) + - U| 
2 2 

(21) 

p(pji Tj) + PjUj = p(j£> 2’ Tg) J°2^2 

(22) 

p l U l = p 2 U 2 

(23) 


Next, energy, momentum, and continuity equations are written for the flow across 
the oblique shock attached to the static wedge whose half -angle is 6. These are 
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(24) 


H(p p Tj) + - (Uj sin 6) 2 = H(p 3 , Tg) + I [U 3 sin(0 - 6)] 2 
2 2 

P (p v Tj) + pj(Uj sin 6 ) 2 = p(p 3 , T 3 ) + p 3 [U 3 sin($ - 5)] 2 (25) 

p 1 U 1 sin 6 = p 3 U 3 sin(0 - 6) (26) 

Since the tangential velocity just upstream of the shock equals the tangential velocity just 
downstream of the shock, there results 

Uj cos e = U 3 cos(0 - 5) (27) 

Next, the density and temperature on the wedge are related to the known pressure on the 
wedge by 


P 3 = P(P S , T 3 ) 


(28) 


The number of equations is reduced by using equation (23) to eliminate U 2 from 
equations (21) and (22) and by using equation (26) to eliminate U 3 from equations (24), 
(25), and (27). In order to more clearly distinguish between the known and unknown 
quantities, a subscripted variable, X R , is used to indicate the unknown quantities. These 
are defined as follows; 


X, = 


I s Pi 


X 5 " p 3 


X 2 = T 1 


X 6 = T 3 


X 0 = 


3 " p 2 


X 4 ” T 2 


X 7 = 


* 

R 


X 8 = 


(Uj sin 9)‘ 
R 


When these substitutions are made, the following eight equations result: 


S(X 3 ,X 4 ) S(p 2 , t , T 2 t ) 

R R 


(29) 


H(X 3 , X 4 ) H(p 2>t , T 2 t ) + 1 X JX ^ 2 q 

R R 


(30) 
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R 

R 

p(x 3 , x 4 ) 

p(x r x 2 ) 

R 

R 

H(x 5 , X 6 ) 

H(Xj, X 2 ) 

R 

R 

P(X 5 , Xg) 

p(x p x 2 ) 



- 1 


l 1 a 7 


(H 


= o 


= 0 



- 1 


L\ X 5 


= 0 


R 


arc sin 




R 


arc tan 


+ X 1 X 3 


I") 


= 0 



x 7 -x 8 


- 6-0 


p(x 5 , x 6 ) _ P3 = 0 


R 


R 


(31) 


(32) 


(33) 


(34) 


(35) 


(36) 


The solution of these equations requires an initial estimate of the X n values that is 
close to the true X n values. The initial values used herein are the result of a perfect- 
gas solution which is discussed as follows: 

Equations (100) and (128) in reference 1 can be combined to yield 


P3 


Ay 


•yy/fy-1) 


P 2,t l(y + 1) 2 J 


(r] - 1 ) 1 ' / '(>'~ 1 )(t; sin 2 6> - 1) 

,y/(y- 1) 


rj' 


where 


(37) 


= 2 y 

y - 1 


Mi 


The following equation is derived from equation (148 (a)) in reference 1 

_2y cos(fl - 6) _ 1_ s ^ n q |- s j n ^ q _ fi) -y sin 6] 

y - 1 f] 2 


(38) 


( 39 ) 


These equations are solved for M-. and 6 . From this, the initial estimates of X ■ 

1 n, 1 
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> 


are determined from equations found in reference 1. These estimates are 


where 



p 2,t 




l/(y - 1 ) 



y + i 


* 5 , 1 ' 


_ y 


— — (M n sin 
- 1 1 


(Mj sin ^ 


ey 




y - l 


x = . 

a 6, i 


^ 2 — (M« sin 6)2 - 1 (M., sin 0) 2 +— - — 

/ - 1 1 J L 1 y - 1 J 




k 2,i 


(m, sin ey 


X 7,i = rM I X 2,i 

x 8,r x v siA 


( 40 ) 


(41) 

(42) 

(43) 

(44) 

(45) 

(46) 

(47) 

(48) 
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(Eqs. (37) to (48) apply only to a perfect gas. The value of y in these equations is 
a perfect-gas value. The values assumed for these calculations are 5/3 for a monatomic 
gas, 7/5 for a diatomic gas, and 4/3 for a gas whose molecules contain three or more 
atoms. ) 

With these values of X n as initial estimates, the procedure for solving equations 
(29) to (36) usually converges rapidly. At this point, flow properties and thermodynamic 
properties can be determined for total and static conditions downstream of the normal 
shock, for static conditions downstream of the oblique shock, and for free -stream condi- 
tions; however, they cannot yet be determined for total conditions in the free stream. 

To do this, the following equations have to be solved: 

s(p ljt . T 1,P = s <Pi* T l> <«> 

H Kt- T i,t> = H K T i> + ^i < 5 °> 


The initial estimate of pj t and Tj t are obtained by assuming the gas to be perfect. 
The perfect -gas relations are 


(Sg-S^/R 

p l,t,i~ p 2,t e 


(51) 


When these initial values are used, the iteration procedure for solving equations (49) and 
(50) converges rapidly. Knowing pj t and Tj t allows the thermodynamic properties 
to be determined for free -stream total conditions. 

Case II - given Pj Tj and p 2 - First, the following equation has to be im- 
plicitly solved for p^ 

P l,t = p ( p l,t’ T 1,P 

Most of the equations that follow are given in case I, but they are repeated here for con- 
venience. The two relations that relate total to static conditions upstream of the shock 
are 


H(p 


s Hf T i,t> = s <Pi, T i> <«> 

l.t- T l.t>- H <Pi, Tjl + iuf (50) 


Next, the energy, momentum, and continuity equations are written for the flow across 
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the normal shock. 


H<Pj, + T 2 ) + iui* 

2 & 

(21) 

p(pj» Tj) + PjUj = P(P2> + p 2^*2 

(22) 

Pl U l = p 2 U 2 

(23) 

Then the relations between total and static conditions downstream of the normal shock are 

S ( p 2,t’ T 2, P “ S ^ p 2’ T 2^ 

(19) 

H ( p 2,t> T 2,P = H ^ p 2, T 2^ + \ U 2 

(20) 

Then the equation that relates the total density and temperature to the known total pres- 
sure downstream of the normal shock is 

p 2, t " p(p 2,t’ T 2,t ) 

(18) 


The number of equations is reduced by using equation (23) to eliminate Ug from 
equations (20) to (22). In order to more clearly distinguish the known and unknown vari- 
ables, a subscripted variable, X Jl , is used to indicate the unknown variables. These are 
defined as follows: 


X l = p l 

X 5 “ p 2, t 

X 2 = T 1 

X 5 = T 2, t 

X 3 = p 2 
X 4 * T 2 

uj 

X 7 = 1 

' R 


When these substitutions are made, the following equations can be written: 


s(x p x 2 ) 




= 0 


R 


R 


( 54 ) 



(55) 


Hftr x 2> T l >t > + 1 


R 


R 


X ? = 0 


h(x 3 , x 4 ) H(X J} x 2 ) j 


R 

+ 

R 

p(x 3 , x 4 ) 

p(x 1? x 2 ) 

R 

R 

S(X 

3 , x 4 ) S(X, 


R 

H(x 3 , X 4 ) 

. H «5' V „ 

R 

R 


+ - X r 



-1 


= 0 


L 1 7 


I-)- 


= o 


R 



= 0 


p< x 5’ X e> oa.t _ A 

R R 


(56) 

(57) 

(58) 

(59) 

(60) 


The solution of these equations requires an initial estimate of the 2^ values that is close 
to the true X n values. A perfect-gas calculation is used to make this initial estimate. 

To do this, the following equation, derived from equation (99) in reference 1, has to be 
solved implicitly for Mji 



Then, when the equations in reference 1 are used, the initial estimate for the X R 
values is 



(62) 


( 63 ) 
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( 64 ) 



(Eqs. (61) to (68) only apply for a perfect gas. ) The use of these values of X n as the 
initial estimates in equations (54) to (60) usually results in the iteration procedure con- 
verging rapidly. From the solution of these equations flow properties and thermodynamic 
properties can be determined for total and static conditions upstream and downstream of 
the normal shock. 


Solution Criteria 

The computer code for these equations includes certain tests to ensure that the so- 
lution is correct. First, all the iterative procedures used to solve the various equations 
and sets of equations have to converge to the prescribed degree. If this convergence is 
not obtained, the solution is flagged; and in some cases, the calculation is terminated. 
Second, even if the solution is mathematically valid, it may not be physically valid. For 
physical validity, the following have to be true: 

(1) All the pressures and temperatures involved in the equations must be within 
range of the state equations. (This also ensures that the fluid is a gas. ) 

(2) The entropy of the gas must increase as the gas passes through a shock. 

If these conditions are not satisfied, the solution is flagged; and in some cases, the cal- 
culation is terminated. 

The following discussion applies only to case I. (This is the case where a wedge - 
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type static -pressure probe is used.) For certain combinations of p 2 t , T 2 p Pg, and 
6, a solutionis obtained which indicates that the oblique shock attached to the wedge is 
the strong shock, rather than the weak shock. Physically, this would be extremely 
unlikely. In fact, the true situation would probably be that the free stream is either sub- 
sonic or that the shock wave is not attached to the wedge. In either case, the solution is 
meaningless. The strong shock is recognized by the fact that a decrease in the wedge 
static pressure is accompanied by a decrease in the free-stream Mach number. Thus, 
the solution is tested by slightly decreasing Pg and solving the equations again; if 
decreases, the oblique shock on the wedge is strong. If this happens, the solution is 
flagged. This test is performed only when the flow downstream of the shock is subsonic. 
Many trial calculations have shown no case in which the downstream real-gas flow was 
supersonic when the shock was strong. 


RESULTS AND DISCUSSION 

For both cases of the preceding analysis, the computer routines are in two sections. 
The routines in the first section include the necessary iteration procedures for solving 
the various sets of equations. The routines in the second section describe the specific 
behavior of the gas. The routines are written in the FORTRAN IV, version 13, language. 

The routines that include the iteration procedures for case I (i. e. , the case where 
the input is p 2 t , T 2 Pg, and 6) are given in appendix B. The routines that include 

the iteration procedures for case II (i.e. , the case where the input is p^ p Tj p and 
P 2 j.) are given in appendix C. The routines for either case have to be used in conjunc- 
tion with one of the sets of routines in appendix D which describe the specific behavior of 
the gas. This appendix includes sets of routines of air, nitrogen, oxygen, normal hydro- 
gen, parahydrogen, helium, argon, steam, methane, and natural gas . These routines can 
also be used in conjunction with the iteration routines in reference 3 to calculate the mass 
flow rate of these gases through nozzles. 

In both cases, the output variables are M^, M^ p , U p Uj p , Mg, and U 2 - The 
second subscript P indicates the value that the variable would have if the gas were per- 
fect. The perfect-gas calculation uses the same input as the real-gas calculation. In- 
volved in this calculation is the perfect-gas specific -heat ratio. The following specific- 
heat ratios were chosen: 5/3 for helium and argon; 7/5 for air, nitrogen, oxygen, nor- 
mal hydrogen, and parahydrogen; and 4/3 for steam, methane, and natural gas. In 
addition, thermodynamic state functions are calculated for static and total conditions up- 
stream and downstream of the normal shock. These are p, p, T, H/R, S/R, C p /R, y, 
and k. Additional variables are calculated for case I. These have to do with the oblique 
shock attached to the static -pressure wedge. These are Mg, Ug, 0, 0 p , and the afore- 
mentioned thermodynamic state functions for the static conditions downstream of the 
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Figure 1. - Ratio of some real-gas flow properties to perfect -gas flow properties for nitrogen - 
case I. Total temperature downstream of normal shock, T^, 300 K; ratio of pressure 
downstream of oblique shock to total pressure downstream of normal shock, p^p^t. 

0.1886; tree-stream Mach number for a perfect gas, /V^ p, 2.5; specific-heat ratio for a 
perfect gas, y p, 1.4; half-angle of static-pressure wedge, 6, 7.5°. 



Figure 2. - Ratio of some real-gas flow properties to perfect -gas flow properties for nitrogen - 
case II. Free-stream total temperature, ^ , 300 K; ratio of free-stream total pressure to 

total pressure downstream of the normal shock, ^ [, 2.004; free-stream Mach num- 
ber for a perfect gas, M} p r 2.5; specific-heat ratio for’a perfect gas, yp, 1.4. 



oblique shock attached to the wedge. 

Figures 1 and 2 are presented to illustrate the deviation of nitrogen from perfect- 
gas behavior. In figure 1, which applies to case I, two ratios are plotted against the 
measured p 2 t - These ratios are Mj/M^ p and The value of Mj p 

is 2.5, and the value of P 3 /P 2 t is constant. For a pressure of 150xl0 5 pascals (new- 
tons per square meter), the Mach number only deviates 0. 8 percent from the perfect- 
gas value, while the mass flux deviates 6. 1 percent. The results are different in fig- 
ure 2. This figure applies to case n. The gas is again nitrogen and the perfect-gas 
Mach number is 2.5. This time pj ^./pg t i s constant. For this case, for a pressure 
of 150x10° pascals, the Mach number deviates 12 percent from the perfect-gas value, 
and the mass flux deviates 5. 2 percent. In both cases, these ratios approach unity as 
the pressure is reduced. This is because the ideal-gas specific -heat ratio of nitrogen 
is sensibly constant at room temperature. This relation is not true for a gas such as 
methane, whose ideal-gas specific heat is variable. Therefore, M^/Mj p and 
PjUj/XpjUjJp would not approach unity, for methane, as the pressure is reduced. 

These routines are designed to be incorporated as a set of subroutines in a main 
computer program. Appendix E includes two typical main programs for the routines 
contained in both appendixes B and C . In addition, the results of sample calculations are 
included. The exact storage requirements depend on the gas. On the average, the rou- 
tines for case I require about 4200 storage locations and the routines for case n require 
about 3500. These stores include a necessary matrix inversion routine. This routine 
is called MINV and is provided by International Business Machines Corporation (IBM) in 
their Scientific Subroutine Package. 


CONCLUDING REMARKS 

Two sets of computer routines are presented that provide a means of calculating the 
properties of a supersonic stream from measured quantities. In the first set, the meas- 
urements are the total pressure and temperature downstream of a normal shock, and 
the pressure on the surface of a wedge -type static -pressure probe. In the second set, 
the measurements are the pressure and temperature in the plenum upstream of a super- 
sonic nozzle and the total pressure behind the normal shock. These routines apply for 
10 different undissociated gases. 

The accuracy of these calculations depends primarily on the accuracy of the state 
equations that describe the gas. For a gas such as air at room temperature and at ple- 
num pressures to 100 bar, the results would probably be accurate to within a percent. 

For a gas such as natural gas, which is a variable mixture of many components, the 
accuracy would not be as good. 

In order to use these routines for gases other than those included in this report, 
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only routines like those included in appendix D would have to be written, 
are straightforward and require only appropriate state equations. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, January 16, 1974, 

501-04. 


These routines 
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APPENDIX A 


SYMBOLS 

C specific heat at constant pressure, J/(kg)(K) 

r 

C v specific heat at constant volume, J/(kg)(K) 

H enthalpy, J/kg 

k isentropic exponent, eq. (16) 

M Mach number 

p pressure, Pa 

R gas constant, J/(kg)(K) 

S entropy, J/(kg)(K) 

T temperature, K 

U velocity, m/sec 

X n unknown quantities in eqs. (29) to (36) where n runs from 1 to 8, and in 

eqs. (54) to (60) where n runs from 1 to 7 

Z compressibility factor 

Zj to Zyj functions of compressibility factor, eqs. (1) to (6) 

a speed of sound, m/sec 

y specific -heat ratio 

5 half -angle of static -pressure wedge 

9 oblique shock angle 

V function of perfect-gas values of and y, eq. (38) 

to £jjj functions of ideal-gas specific heat, eqs. (7) to (9) 

3 

p density, kg/m 

Subscripts: 

i initial estimate of a variable in an iteration process 

id ideal-gas condition 

P perfect-gas condition 

t total conditions 
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1 refers to free -stream conditions except when appended to X 

2 refers to conditions downstream of normal shock except when appended to X 

3 refers to conditions downstream of oblique shock except when appended to X 
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APPENDIX B 


DESCRIPTION AND CARD LISTING OF COMPUTER ROUTINES 
THAT APPLY TO CASE I 

Case I applies to local measurements of stagnation pressure Pg stagnation tem- 
perature Tg and static pressure p^ on a wedge of half -angle 6. 

In order to make a workable set of routines, the routines in this section have to be 
combined with one of the sets of routines in appendix D. The particular set that is cho- 
sen depends on the type of gas. This combined set of routines is referenced in the main 
program by the following statement: 

CALL RWEDG(PA2, TA2, PB3, DELTAS) 

All the output is returned through labeled COMMON. Therefore, this COMMON 
statement has to be in the main program: 

COMMON/OUTW/TP(8, 5), VM(10), CONV(12), Z(6, 5),KODE(ll) 

There are certain variables that have to do with the convergence requirements of 
the various iteration procedures in these routines. These variables are initialized in a 
BLOCK DATA subprogram. However, if the following COMMON statement is put in the 
main program, the value of these variables can be altered: 

COMMON/LIMITW/EDA2, EMAT, EDTAl, EPWD 

The symbols that apply to these routines are defined as follows: 

PA2 Total pressure downstream of normal shock, pg Pa 

TA2 Total temperature downstream of normal shock, Tg ^ K 

PB3 Pressure downstream of oblique shock attached to wedge -type static - 

pressure probe, pg, Pa 

DELTAS Half-angle of wedge-type static -pressure probe, 6, deg 

The array TP(I,J) contains thermodynamic state functions. The subscript I iden- 
tifies the function, and the subscript J identifies the region. The function definitions 
for the first subscript are as follows: 

TP(1,J) Pressure, p, Pa 

TP(2,J) Density, p, kg/m^ 
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TP(3,J) Temperature, T, K 

TP(4,J) Enthalpy, H/R, K 

TP(5,J) Entropy, S/R 

TP(6, J) Specific heat, Cp/R 

TP(7,J) Spe cific -heat ratio, y 

TP (8, J) Isentropic exponent, k 

The region definitions for the second subscript are as follows: 

TP(1, 1) Total conditions in free stream 

TP(I, 2) Free-stream conditions 

TP(I, 3) Total conditions downstream of normal shock 
TP(I,4) Conditions downstream of normal shock 

TP (I, 5) Conditions downstream of oblique shock attached to wedge -type static - 
pressure probe 

(Thus, TP(3, 5) would be the temperature downstream of the oblique shock and TP(8, 2) 
would be the free-stream isentropic exponent.) 

VM(1) Free-stream Mach number, Mj 

VM(2) Perfect-gas calculation of free-stream Mach number, M^ p 

VM(3) Free-stream velocity, Uj, m/sec 

VM(4) Perfect-gas calculation of free-stream velocity, U-, m/sec 

VM(5) Mach number downstream of normal shock, Mg 

VM(6) Velocity downstream of normal shock, Ug> m/sec 

VM(7) Angle of oblique shock attached to wedge, 8 , deg 

VM(8) Perfect-gas calculation of oblique shock angle, 0 p , deg 

VM(9) Mach number downstream of oblique shock, Mg 

VM(10) Velocity downstream of oblique shock, Ug 

CONV(l) to CONV(8) 
where 

X - X 

CONV(l) = -iii Li ' 1 

x t „ 1 
I, n-l 
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C ONV (9) to CONV(12) 
where 


CONV(9) = 1 


2,t 


Z 2,t p 2,t RT 2,t 


CONV(IO) = 


CONV(ll) = 


( T i t ) - (Tj J 


(T 1 i 

i > 1 n -1 


CONV(12) = 



The array Z(I> J) contains the compressibility factor functions as represented by 
equations (1) to (6). The subscript I identifies the function, and the subscript J iden- 
tifies the region. The second subscript definitions are the same as those for the array 
TP. The first subscript definition is 


Z(1,J) . . . Z (6, J) for Zj to Z yI 


The following symbols represent integers to indicate various error conditions. If 
all the integers are zero, a valid calculation has been performed. If the integers are 1, 
errors exist. These errors are 


KODE(l) 
KODE (2) 

KODE(3) 

KODE (4) 


Equals 1 if p^ t or Tj t is out of range of the state equations. 

Equals 1 if pg ^ or T 2 ^ is out of range of the state equations. A value 
of 1 terminates the calculation. 

Equals 1 if pj or Tj is out of range of the state equations. A value of 1 
terminates the calculation. 

Equals 1 if pg or Tg is out of. range of the state equations. A value of 1 
terminates the calculation. 
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KODE (5) 

KODE(6) 

KODE (7) 
KODE (8) 

KODE (9) 

KODE (10) 
KODE (11) 
EDA 2 


EMAT 


EDTA1 

EPWD 


Equals 1 if pg or Tg is out of range of the state equations. A value of 1 
terminates the calculation. 

Equals 1 if the iterative procedures for the perfect-gas calculations for 

M 1 P’ M 2 P’ and 9 P fail to conver £ e - 
Equals 1 if the iterative procedures for calculating P 2 ^ fails to converge. 

Equals 1 if the iterative procedures for calculating Xj to Xg fail to con- 
verge . 

Equals 1 if the iterative procedures for calculating p^ ^ and ^ fail to 
converge . 

Equals 1 if Sj > S£ or if > Sg. 

Equals 1 if the shock attached to the wedge is the strong shock. 

Maximum value of |CONV(9)| permitted. This is initialized at a value of 
lxlO -6 

8 

Maximum value of V' JCONV(I)| permitted. This is initialized at a 

1=1 

value of lxicr^. 

Maximum value of |CONV(10)| + |CONV(ll)| permitted. This is initialized 
at a value of 2xl0~ 6 . 

Maximum value of |CONV(12) | permitted. This is initialized at a value of 
1X10' 6 . 


The routines that apply to case I are described briefly here. These routines are 
identified by their deck names. 


Deck RGWED 

This is the subroutine in which the output variables TP and VM are calculated. 
The iteration procedures for solving equations (18), (29) to (36), and (49) and (50) are in 
this routine. 


Deck RGDFW 

In this subroutine, the inverse of the derivative matrix, which is required to solve 
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equations (29) to (36), is calculated. In addition, this routine requires the matrix inver- 
sion subroutine MINV supplied by IBM in their Scientific Subroutine Package. 


Deck PGWED 

This subroutine is used to calculate Mj p, M 2 p, and 0p. The iterative proce- 
dures necessary to solve equations (37) to (39) are in this subroutine. 


Deck WDATA 


This BLOCK DATA subprogram contains the initial values of the convergence limits 
for the various iterative procedures. 

A workable set of routines includes not only these decks, but also one of the sets of 

decks in appendix D. The card listing of these decks follows. 

C DECK RGWEO 
C 

C THIS IS THE SUBROUTINE WHERE THE FLOW FUNCTIONS AND THERMODYNAMIC 
C STATE FUNCTIONS OF A SUPERSONIC STREAM ARE CALCULATED FROM MEASURED 
C QUANTITIES. THESE QUANTITIES ARE THE TOTAL PRESSURE AND TOTAL 
C TEMPERATURE DOWNSTREAM OF A NORMAL SHOCK, ANO THE PRESSURE INDICATED 
C BY A WEDGE STATIC PRESSURE PROBE. 

C 

SUBROUTINE RWEDG ( PA2 , TA2, PB3, DELTAS ) 

COMMON /OUTW/ TP < 8, 5 ) , VM( 10 } ,CONVU 2 ) , Z ( 6, 5 > , KODE ( 11 ) 

COMMON /ROW/ W(8,8),X(8) ,X13, SQX13,X15, SQXl 5, V7, W7, CVB1 ,CVB2, CVB3 
COMMON /GDATW/ G1,G2 , G3, G4, G5, G6, G7, G8, G9, G10 
COMMON /LDATA/ XKV, R , XMW ,RC , 02 , G 
COMMON /LIMITW/ E0A2 , EMAT, EDTA1 , EPWD 

DIMENSION ZA1(6), ZB1(6), ZA2<6), ZB2(6J, ZB3(6), F<8>, 0UT(1D3>, 
1SAVE ( 1 33 ) 

EQUIVALENCE ( TP( 1) , 0UT( 1 )) , < ZAM 1 > ,Z < 1 )) » IZBim,Z(7>), (ZA2C1), 

1Z ( 13 ) ) * ( ZB2 ( 1 ) , Z ( 19 ) ) , { ZB3 1 1 1 , Z ( 25 > ) 

DOUBLE PRECISION CP, CS, CH, DS A2, DHA2, DSA1, DHA1, DSB1.DHB1 ,DS82,DHB2 , 
1DHB3 

LOGICAL LGFN 

DATA RADI AN/. 0174532925/ 

00 1 1=1,103 

SAvem-o.o 

1 0UT(I>=0.0 
DO 2 J = 1 , 5 
00 2 1=1,3 

2 Z(I,J)=1.0 
DO 3 1=1,8 

xm=o.o 

3 F ( I ) =0.0 
TP(1,3»=PA2 
TP(3,3)=TA2 
T P ( 1 , 5 ) =PB3 
NSAVE=D 
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THE PERFECT GAS VALUES OF THE SHOCK ANGLE* AND THE MACH NUMBER, BOTH 
UPSTREAM AND DOWNSTREAM OF THE NORMAL SHOCK ARE CALCULATED IN PWEDG. 

CALL PWEDG (TP( 1, 5 ) /TP< 1 , 3 ) , DELTAS, VM( 2 ) , XM2, VM< 8 ) ,CONV< 12) ,K0DE(6 

1 ) > 

IF (K00EI6) .EQ.l ) RETURN 
DELT A=DELTAS*RADI AN 
SQM1=VM( 2 )**2 

SQSM1=SQM1*SIN(VM(8)*RADIAN) **2 

THE ITERATION PROCEDURE FOR CALCULATING THE TOTAL DENSITY DOWNSTREAM 
OF THE NORMAL SHOCK STARTS HERE. 

A = TPU,3>/<R*TPI3,3) ) 

TP( 2* 3 ) =A 

IF (LGFNITPI 1,3) ,TPt 3,3) ,KODE( 2 ) , 2 A2 ) ) RETURN 
DO 5 1=1,20 

CALL ZETA ( 1,TPI2,3) ,TP( 3,3) ,ZA2) 

CONV<9)=1.0-TP{ 1,3)/ (ZA2(1 ) *R*TP< 2, 3 1 *TP{ 3, 3) ) 

IF ( ABS ( CONV (9) ).LT.EDA2) GO TO 6 

TP( 2 ,3 )*TP( 2 ,3)*( 1.0-1 ZA2( 1 ) -A/ TP I 2»3))/ZA2(3) ) 

KODE 1 7 ) = 1 
RETURN 

CALL ZETA 13,TP12,3) ,TP( 3,3) ,ZA2> 

IF 1LGFN1TP1 1,3),TPI 3,3) ,KODE( 2), ZA2) ) RETURN 
DSA2=CS(TP(3,3) )-DLOG( DBLE l TPI 2, 3) ) )-ZA2(4) 

TP(5,3)=DSA2 

DHA2=CH1TP(3»3))+TP13»3)*<ZA211)-ZA2(5>) 

TP(4,3)=DHA2 

CV=CP(TP13,3))-ZA2(6) 

GA=ZA2 (3)+ZA2t2)**2/CV 
TPl 7, 3 )=GA/ZA2.1 3) 

TP(8,3)=GAZZA2(1) 

TP(6,3)=CV*TP17,3> 

S0M2=( SQH1+G5)/(G7*SQM1-1.0) 

THE INITIAL ESTIMATES OF THE VALUES OF X ARE CALCULATED BY THE 
FOLLOWING STATEMENTS. 

A=1.0*Gi*SQM2 
X(3)=TP(2»3)/A**G3 
X ( 4 ) =TP( 3 , 3 ) / A 

X ( 1 ) =X 1 3)*( SQM1+G5 ) / ( G6*SQM1 ) 

X(2)=TP(3,3)/I1 . 0+Gl*SQMl) 

A=G6*$QSMl/ < SQSM1+G5 ) 

X (5 ) = A*X( 1 ) 

X(6)=X12>*IG7*SQSM1-1.0>/1G6*A) 

X 1 7 ) =G2*SQM1*X ( 2 ) 

X I B ) =G2*S0SM1*X(2) 

VMC 4 ) = SQRT 1 R*X(7 ) ) 

THE ITERATION PROCEDURE FOR CALCULATING THE VALUES OF X STARTS HERE. 
DO 11 1=1,20 

CALL ZETA { 3, X (1) ,X1 2 ) , ZB1 ) 

CALL ZETA ( 3 , X ( 3 ) , X ( 4 ) , Z82 ) 



a> -si 


CALL ZETA ( 3 , X ( 5 > ,X< 6 J , ZB3 > 

TP(i,2) = ZBl(l)*Xm*X(2)*R 
TP(1,4)=ZB2(1>*X(3>*X<4)*R 
IF (LGFN(TP(1,2),X<2),K0DE(3),ZB1) ) RETURN 
IF ( LGFNI TPI 1,4) »X( 4 ) ,KODEI 4 ) ,ZB2 ) ) RETURN 
IF <LSFNITP<l,5),X<6),K0DE<5>,Z83n RETURN 
DSB2=CSIX(4) >-DL0G(0BLE(X(3> ))-ZB2(4) 
DHB1=CHI X(2) )*X(2)*{ZB1(1)-ZB1(5) ) 
DHB2=CM(X(^))+X<4)*(ZB2( 1J-ZB2I5) ) 

DHB3-CHI X ( 6 ) KX(6)*I ZB3I 1)-Z83<5> > 
xi3=xm/xm 
SQX13=X13**2 
X15»XIU/X(5) 

SQX15*X15**2 

V7=SQRT(SQX15*X(8)/(X(7)-X(8)) ) 
W7=SQRTIXI8>/X(7)) 

Fill =0SB2-DS A2 

F ( 2 ) =DHB2-0HA2+SQX13*X (71/2.0 
F<3)=DHB2-DHB1+X<7>*< SQX13- 1 .0 ) /2. 0 
F (4)=DHB3-DHBl>X(8)*(SQXI5-1.0)/2.0 
Fl=TP(l,2)/R 
F3=TP(l,4)/R 
F5-XI5)*X(6J*ZB3(i> 

F (51 =F3-F1+X(1)*X17) *1X13-1*0) 

F <6 )=F5-F1+X(1)*X( 81*1X15-1.0) 
F(7)=ARSIN(W7)-ATANI V7) -DELTA 
F(8)=F5-TP< I,5I/R 
IF ( MODI 1,5 ) «£Q. I ) CALL RDFW 
DO 8 J=1 , 8 
CONV(J)=0.3 
DO 7 K*l,8 

CONV C J ) =CONV (J)+W(J,K)*F(K) 

CONV I J ) =CONV t J ) /X( J) 

TEST=0.0 
DO 9 J-l, 8 
9 TEST=TEST+ABS(CONVI J )) 

IF { TEST. LT. EMAT ) GO TO 12 
DO 10 J=1 ,8 

id xi j)=xij)*( i.o+conv( j n 

11 continue 

c 

KOBE C 8 I = 1 

12 CVB1=CP(X(2))-ZB1I6) 
GA=ZB1I3)+ZB1I2)**2/CVB1 
TP( 7 ,2 ) *GA/ZBl( 3 ) 

TP( 8,2)=GA/ZB1( 1 ) 

TP(6,2)=*CVB1*TP< 7,2) 

TPI 3,2 ) *X (2 ) 

TP(2,2)=X(l) 

TP( 4,2)=DH81 

TPI5,5)=CSIX(6) )-ALOG(XI 5) >-ZB3<4) 

TP( 4,5 )=DHB3 
TP<2,5)=X15) 

TPI3,5)=X{6> 

GA*ZB3(3)+ZB3C2)**2/CVB3 
TPI 7 ,5 )=GA/ZB3 ( 3 ) 

TP ( 8 , 5 ) *GA/ ZB3 1 1 ) 

TP(6,5)=CVB3*TPI7,5) 

VM(3)^S0RT(R*X(7) ) 
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VM( l)=tfM{3) /SORT ( TP ( 8,2>*TP( 1,2)/TP(2,2H 
VM(7)=ARSIN(SQRT(X(B)/X(7) ) ) 

VMIIO) = COS( VM ( 7 > )*SDRT(R*X(7) ) /COS ( VM ( 7 ) -DEL TA > 
VM(9)=YM(1D)/SQRT(TP(8,5)*TP(1,5)/TP<2,5) > 

IF (NSAVE.NE.l) GO TO 14 
IF(VM(1).LT.SAVE(41) JKODEI 1 1 J = 1 
DO 13 1=1,102 

13 OUT ( I ) =SAVE ( I ) 

RETURN 

14 D$B1=CS(TP(3,2> )-DL0G( D8LE< TP( 2, 2) ) J-ZBK4) 

TP(5,2 MDSB1 
VM(7)=VM<7) /RADIAN 
TP ( 3 ,4 ) =X ( 4 ) 

TP ( 2 ,4 ) =X ( 3 ) 

TP ( 4,4 ) =DHB2 
TP< 5,4)=DSB2 
CVB2=CP(TP( 3,4) ) -ZB2 C 6 ) 

GA=Z82(3)+ZB2(2)**2/CVB2 
TP(7,4)=GA/ZB2<3) 

TP ( 8 , 4 ) =GA/ Z82 ( 1 ) 

TP(6,4)=TP(7,4) *CVB2 
VM(6)=SQRT(2.0*R*(TP(4,3)-TP(4,4) ) ) 
VM(5)=VM(6)/SQRT(TP<8»4)*TP(1*4)/TP(2,4) ) 

IF <TP(5,2).GT.TP(5,4).OR.TP(5,2).GT.TP<5,5> > KODE(10)=1 

THE ITERATION PROCEDURE FOR CALCULATING THE TOTAL DENSITY AND 
TEMPERATURE STARTS HERE. 

TP(3,1)=TP(3,3) 

TP (2,1 )=TP(2,3)*EXP(SNGL(DSB2-DSB1) ) 

DO 15 1=1,20 

CALL ZETA ( 3 , TP ( 2, 1 ) , TP( 3, 1 ) , Z All 
TP(1,1)=TP(2,1)»TP(3,1)*R*ZA1<1) 

IF (LGFN(TP(i,l) ,TP< 3,1) ,K0DE(1),ZA1) ) K00E<1)=1 
DSA1=CS(TP(3,1))-DL0G(DBLE(TP(2, 1) ) )-ZAl(4) 

F 1=DSB1— DSA1 

DHAl=CH(TP(3,l) )+TP(3,l)*(ZAl(l)-ZAl(5) ) 

F2=X(7 ) /2 . D+DHB1-DHA 1 
Fll=-ZAi(2)/TP(2,l) 

CVA1=CPCTP(3,1) )-ZAl(6) 

F12=CVA1/TP(3,1) 

F21=TP(3,1)*(ZA1(3)-ZA1( 2) )/TP(2,l ) 

F22=CVA1+ZA1<2) 

DET = FU*F22-F12*F21 

C0NV(1D) = (F1*F22*-F2*F12)/(TP(2,1)*0ET) 1 

CONV< 11) =<F2*Fll-Fl*F21)/( TP( 3, 1 )*OET> 

TEST = ABS(CONV(lO> M-ABSiCONVUll) 

IF ( TEST.LT.EDTA1) GO TO 16 
TP (2,1 )=TP( 2 , 1 ) *( 1. 0+CONV( 10) ) 

15 TP(3,1)=TP(3,1)*(1.0+CONV(11>) 

C 

KODE ( 9 ) = 1 

16 TP ( 4, 1 )=DHA1 
TP ( 5 , l ) =DS A1 

GA=ZA1(3)+ZA1(2)**2/CVA1 

TP(7,1)=GA/ZA1(3) 

TP(8,1 )=GA/ZA1 (1 ) 

TP(6,1)=CVA1*TP(7,1) 

IF ( VM(9) .GT.l.D) RETURN 
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NSAVE=1 
DO 17 1=1,103 
17 SAVe(I)=OUT(I) 

TP( 1 ,5>=TPI 1*5 >-0.001*(TP( 1,3)-TP(1,5)) 

GO TO A 

END 


OECK RGDFW 

THE INVERSE OF THE DERIVATIVE MATRIX REQUIRED FOR DETERMINING THE 
VALUES OF X IS CALCULATED IN THIS SUBROUTINE, 

SUBROUTINE RDFW 

COMMON /ROW/ W<8,8), X(8> ,X13,SQX13,X15,SQX15,V7, W7,CVB1,CVB2,CVB3 
COMMON /OUTW/ E < 62 ) , Z A1 ( 6) , Z B1 ( 6 ) , Z A2 ( 6 ) , ZB2 ( 6 ) ♦ Z 83 1 6 ) , KODE ( 1 1 I 
DIMENSION XW<6)> YW<8> 

DOUBLE PRECISION CP 
DO 1 1=1,8 
DO 1 J=l, 8 
1 W<I,J)=0.0 

SV7=1.0+V7**2 
SW7=SQRT( 1.0-W7**2> 

CVB1 =CP(X ( 2 ) 1-ZB1 < 61 
CVB2=CP<X<A) )-ZB2(6) 

CVB3=CP(X<6» I-Z83I6) 

W(1,3)=ZB2(2)/X(3) 

W { l , A) =-CVB2/X( A ) 

W (2 , 1 ) =— X 1 7 ) *SQX13/X ( 1 ) 

W ( 2 , 3 ) = (SQX13*X< 7)-X (AJ*(ZB2I3I-ZB2(2>IJ/X(3) 

W(2,A)=-CVB2-ZB2<2) 

W(2,7)=-S 0X13/2,0 
A=X(2)*(ZB1(3)-ZB1(2 ) ) 

W(3,n=(A-SQX13*X(7> >/X(l) 

W(3,2)=CVB1+ZB1C2) 

W(3,3)=W(2»3) 

W ( 3, A) =W( 2, A ) 

W(3»7)=0.5+W(2,7) 

WIA, l I *( A—XI8 )*SQX15 )/X( 1) 

WIA, 2)=VM3,2) 

W(A,5)=(SQX15*X(8>-X(6)*(ZB3(3)-ZB3(2) H/XI5) 

W(A»6)=-CVB3-Z83<2> 

W(A,Bl=<l,0-SQX15)/2.0 

Wt 5 » 1 ) *X( 2 1 *ZBH 3>-X 1 7 )*( 2 ,0*X13-1, 0) 

W(5,2)=X< 1)*ZB1(2) 

W(5,3)«XI7)*SQX13-X( A»AZB2(3J 
W ( 5, A) *-X ( 3 ) *ZB2 ( 2 > 

W(5,7)=U.0-X13)*Xm 

W(6,1)=X(2I*ZB1I3)-X( 8)*(2,0*X15-1,0) 

W (6 , 2 ) =W( 5,2) 

W<6,5) *X(8>*SQX15-X< 6>*ZB3( 3 ) 

W(6,6)=-X<5)*ZB3(2) 

W(6,8)=Xm*(1.0-X15) 

W(7,l)=V7/( SV7AXC1)) 

W(7»5)=-WI7,1)*X15 
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W<7,7)=W7/t 2.0*X<7)*SW7)-V7/(2.0*SV7*<XI T J -X { 8 » ) ) 
WI7,8)=(V7/(W7*SV7*(Xm-X<8)l)-l.0/IXI7)*SW7> )/< ?.0*W7> 
W(8»5)=-X(6)*ZB3(3) 

W(8»6)=-X(5)*ZB3t2) 

CALL MINV (W,8,DET,XW,YW ) 

MINV IS THE SUBROUTINE USED TO INVERT THE DERIVATIVE MATRIX. THIS 
SUBROUTINE IS PART OF THE SCIENTIFIC SUBROUTINE PACKAGE SUPPLIED 
BY I.B.M. 

RETURN 

END 


C DECK PGWED 
C 

C THE PERFECT GAS VALUES FOR VELOCITY, SHOCK ANGLE, AND THE MACH NUMBER 
C BOTH UPSTREAM AND DOWNSTREAM OF THE NORMAL SHOCK ARE CALCULATED IN 
C THIS SUBROUTINE. 

C 

SUBROUTINE PWEDG ( RAT 1 0, DELTA, XM1, XM2, THETA, DELR , KK ) 

COMMON /L DAT A/ AI5),G 

COMMON /GDATW/ G1 , G2 , G3, G4 , G5, G6, G7, G8, G9, G10 
COMMON /LIMITW/ EDA2 , EM AT, EDTA1 , EPWO 
DIMENSION DEL ( 3 I 

DATA RADI AN, PI, KG/. 01 74532925, 1.5707963,0/ 

FF ( A ) = ( 1. 0-1 .0/ A >**G3* (ST2-1.O/AJ/GI0-R 
XX ( A)=G8*C0SIA-D)/(ST*<SIN< 2.0*A-D ) -GSD ) > 

FFMI A)=SIN( A)**2/G10-R 
IF (KG.NE.O) GO TO 1 
Gl-IG-l.O >/2.0 
G2 = G 

G3*1.0/(G-1.0) 

G4=G*G3 

G5=2.0*G3 

G6= C G+l .0 )*G3 

G7=2.0*G4 

G8=2.0*G7 

G9=G6**G6 

G10=I(G+1.0)**2/I4.0*G) ) **G4 
KG= 1 

1 KKK=0 

DR=O.0 
R=RATI 0 

D=OELTA*RADI AN 
GSD=G2*SIN< D> 

IF (D.LE.O.O.OR.GSO.GT.l.O) GO TO 11 

TMI N=( D+ARS l N( GSD) )/2.0 

TMAX=PH-D-TMIN 

Tl=TMlN+l«0E-4 

F1=FFMI TMIN ) 

IF ( FI • GT. 0 . O.OR.FFM (TMAX) .LT.O.O) GO TO 11 
DELI 1) s(TMAX-TMIN-2.0E-4)/35.0 
DELI2) =5.0*DEL(1> 

OEL I 3) =DEL( 1 > 

DO 2 J=l,3 
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DO 2 L =1 * 5 
T2=T1*DEL( J) 

ST= S IN(T2 ) 

ST2 = $T **Z 
X=XX C T 2 ) 

F2 = FF< X) 

IF (ABS(F1*F2)»LT«ABS(F1)+ABS(F2) ) SO TO 3 
F I = F2 
T1 = T2 
GO TO 11 

3 I N=0 

00 9 1=1,23 
DR=F2/R 

IF lABS(DR).LT.EPWD) GO TO 10 
DT* ( T2-T1 »/(Fl/F2-1.0) 

T1 = T2 
F1 = F2 
T2=T2*DT 

IF < T2-TMIN-1.0E-4) 4,4,5 

4 T2=TMIN+1.0E-4 
GO TO 7 

IF (TMAX-T2-1.0E-4) 6,6,8 
T2=TMAX— l.OE-4 
IF < IN.EQ.l ) GO TO 11 

1 N= 1 

B ST=S IN( T2 ) 

ST2 = ST **2 
X = XX < T2 ) 

9 F2 = FF( X) 

KKK = 2 

10 IF <X.GT.G7> GO TO 12 

11 KKK=1 
X=0.0 
T2=0. 0 
X2=0. 0 

12 OELR=DR 

IF (X.EQ.O.OI GO TO 13 
SX=ST2*X 

X2=(SX/G7+G5)/(SIN(T2-D)**2*(SX-1«0) ) 

13 KK=KKK 
THETA=T2/RADIAN 
XM1=SQRT f X/G7 ) 

XM2 = SQRT { X2 I 
RETURN 

END 


C DECK WDATA 


BLOCK DATA 

COMMON /LIMITW/ E(4) 

DATA E/1.0E-6,1.0E-5,2.0E-6, 1.0E-6/ 
END 
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APPENDIX C 


DESCRIPTION AND CARD LISTING OF COMPUTER ROUTINES 
THAT APPLY TO CASE E 

Case II applies to measurements of nozzle flow using the upstream total pressure 
Pj t , the upstream total temperature Tj t , and the stagnation pressure pg t at the 
nozzle exit. 

In order to get a workable set of routines, the routines in this appendix have to be 
combined with one of the sets of routines in appendix D. The particular set that is 
chosen depends on the type of gas. This combined set of routines is referenced in the 
main program by the following statement: 

CALL RNS(PA1, TA1, PA2) 

All the output is returned through labeled COMMON. Therefore, this COMMON 
statement has to be in the main program: 

COMMON/OUTNS/TP(8,4), VM(6),CONV(9),Z(6,4),KODE(8) 

There are certain variables that have to do with the convergence requirements of 
the various iteration procedures in these routines. These variables are initialized in a 
BLOCK DATA subprogram. However, if the following COMMON statement is put in the 
main program, the value of these variables can be altered: 

COMMON/LIMITN/EDA1, EMAT, EPNS 

The symbols that apply to these routines are defined as follows: 

PA1 Free -stream total pressure, p^ t , Pa 

TAl Free -stream total temperature, Tj ^ K 

PA2 Total pressure downstream of a normal shock, pg ^ , Pa 

TP(I, J) Except for the fact that J does not assume the value of 5, the definition of 

the elements in this array is the same as that for the TP array in 
appendix B. 

VM(I) Except for the fact that I does not exceed 6, the definition of the elements 

in this array is the same as that for the VM array in appendix B. 

CONV(l) to CONV(7) 
where 
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C0NV(I) = 


X I, n ~ X I, n-1 
X I, n-1 


CONV(8) and CONV(9) 
where 


Z(I,J) 



Except for the fact that J does not assume the value of 5, the definition of 
the elements in this array are the same as that for Z in appendix B. 


The following symbols represent integers to indicate various error conditions. If 
all the integers are zero, a valid calculation has been performed. If the integers are 
unity, errors exist. These errors are 


KODE(l) 

KODE(2) 

KODE (3) 

KODE(4) 

KODE(5) 
KODE (6) 
KODE (7) 

KODE (8) 
EDA1 


Equals 1 if pj t or Tj t are out of range of the state equations. A value 
of 1 terminates the calculation. 

Equals 1 if p 2 t or Tg t are out of range of the state equations. A value 
of 1 terminates the calculation. 

Equals 1 if pj or T^ are out of range of the state equation. A value of 1 
terminates the calculation. 

Equals 1 if p 2 or T 2 are out of range of the state equation. A value of 1 
terminates the calculation. 

Equals 1 if the iterative procedures for calculating Mj p fail to converge. 

Equals 1 if the iterative procedures for calculating pj t fails to converge. 

Equals 1 if the iterative procedures for calculating Xj to Xy fails to 

converge. 

Equals 1 if > S 2 

Maximum value of |CONV(8)| permitted. This is initialized at lxl(T 6 . 
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EMAT 


7 

Maximum value of 

1=1 



| C0NV(I) I . 


This is initialized at 1x10 


-5 


-6 

EPNS Maximum value of |CONV(9) j permitted. This is initialized at 1x10" . 

The routines that apply to case n are described briefly here . These routines are 
identified by their deck names. 


Deck RGNS 

This is the subroutine in which the output variables TP and VM are calculated. 
The iteration procedure for solving equations (53) and (54) to (60) are in this routine. 


Deck RGDFNS 

In this subroutine, the inverse of the derivative matrix is calculated. This is re- 
quired to solve equations (54) to (60). In addition, this routine requires the matrix in- 
version routine MINV supplied by IBM in their Scientific Subroutine Package. 


Deck PGNS 

This subroutine is used to calculate M^ p- The iterative routine used to solve equa- 
tion (61) is in this subroutine. 


Deck NSDATA 

The BLOCK DATA subprogram contains the initial values of the convergence limits 
for the various iterative procedures. 

A workable set of routines has to include not only these decks, but also one of the 
sets of decks in appendix D. The card listing of these decks follows. 
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C DECK RGNS 
C 

C THIS IS THE SUBROUTINE WHERE THE FLOW PROPERTIES A 'ID THERMODYNAMIC 
C STATE FUNCTIONS OF A SUPERSONIC STREAM ARE CALCULATED FROM MEASURED 
C QUANTITIES. THESE QUANTITIES ARE THE PRESSURE AND TEMPERATURE IN A 
C PLENUM UPSTREAM OF A SUPERSONIC NOZZLE, AND THE TOTAL PRESSURE 
C DOWNSTREAM OF A NORMAL SHOCK AT THE NOZZLE EXIT. 

C 

SUBROUTINE RNS ( PA1, TAI , PA2 > 

COMMON /OUTNS/ TP ( B, 4 1 , VM< 6 ) ,CONV I 9 > , Z< 6, A > , KODE ( 8 > 

COMMON /RONS/ W ( 7 , 7) , XI 7 ) , X13, SQX 1 3 
COMMON /GDATNS/ G1,G2,G3,G4,G5,G6,G7,G8, G9,G10,G11 
COMMON /LDATA/ XKV, R , XMW ,RC, D2, G 
COMMON /LIMITN/ EDA1 , EMAT , EPNS 

DIMENSION IAII6I, ZB 1 1 6 ) , ZA2I6), ZB2I6), F(7), 0UTI79) 

EQUIVALENCE I TPI 1 } , OUT (II), (lAllll.ZdH, (ZBim»Z(7M, (ZA2(I), 
1 Z 1 1 3 ) } , IZB2<1>,Z(19>> 

DOUBLE PRECISION CP, CS,CH, DS A1 ,DSA2, DSBl , DSB2 ,DHA 1 ,OHA2, DHB1 ,DHB2 
LOGICAL LGFN 
DO 1 1*1,79 

1 OUT (I) *0.0 
DO 2 J = l,4 
DO 2 1*1,3 

2 ZIl,J)*1.0 
DO 3 1=1,7 
XII )*0.0 

3 Ftl>*0.0 
TP(1,1I*PA1 
TP ( 3 , i ) =T A1 
TP 1 1 , 3 ) =P A2 

THE PERFECT GAS VALUES OF THE MACH NUMBER, BOTH UPSTREAM AND 
DOWNSTREAM OF THE NORMAL SHOCK ARE CALCULATED IN PNS. 

CALL PNS (TP ( l , 3 )/TP 1 1, 1 ) , VMi 2 ) »VM( 5 ) ,CONV( 9 ) ,KODE I 5 ) } 

IF ( LGFNI TP(1,1),TP(3,1I , KODEI 1 ), Z A1 ) ) RETURN 

THE ITERATION PROCEDURE FOR CALCULATING THE PLENUM DENSITY STARTS 
HERE. 

A*TP(i,l)/| R*TPI 3,1) ) 

TP ( 2 , 1 ) =A 
DO A 1=1,20 

CALL ZETA ( 1,TP ( 2, 1 I , TPI 3, i ) ,ZA1) 
CONV(8l=1.0-TP(l,l)/CZAim*R*TP(2,l»*TP(3,in 
IF ( ABSI CONV (8) I.LT.EDA1 ) GO TO 5 
AA=(ZA1(1)-A/TP(2,1> I/ZAM3) 

TP (2 ,1 ) “TPI 2 » 1 )*( l.D-AA) 

KODE ( 6 ) =1 

5 CALL ZETA I 3, TPI 2, 1) , TPI 3, 1) , ZA1 ) 

IF (LGFNI TP( 1, 1 ) , TPI 3, 1 ) ,KODE( 1 ) > ZA1 > } RETURN 
DHA1*CHITPI 3,1) ) *-TP( 3, 1 > *1 Z A1 1 1 ) -Z A1 1 5 )) 

TPI A, l ) =DHA1 

DSA1*CS(TPI3,1) ) -DLOGI OBLEI TPI 2, 1 1 ) J-ZAKA) 

TPI5,U=DSA1 

CV=$NGL(CP(TP(3,i)))-ZAl(6) 

GA=ZA1(3)+ZA1(2) **2/CV 
TPI7,1»*GA/ZA1I3) 

-TP < 8 , 1 )=GA/ZA1 1 1 ) 
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TP(6,1)=CV*TPI7,1) 

IF ( KODE I 5 ) . EQ. 1 1 RETURN 

THE INITIAL ESTIMATES OF THE VALUES OF X ARE CALCULATED BY THE 
FOLLOWING STATEMENTS. 

SQM1=VM<2 >**2 
AA=(1.0+G1*SQM1I 
X { 1 ) =A/AA**G3 
X ( 2 ) =TP( 3 » 1 ) / AA 
AA=G6*SQM1/ (SQMH-G5) 

X(BI=X(U*AA 
AAA=(G7*SQM1-1.0)/G6 
X(4)=X(2)*AAA/AA 
X< 5)=A*AA**G4/AAA**G3 
X(6I=TP(B»l) 

X 17 ) =G2*X < 2 ) *SQMI 

VM<4l=SQRT(Xm*R) 

THE ITERATION PROCEDURE FOR CALCULATING THE VALUES OF X STARTS HERE. 

00 10 1 = 1,20 

CALL ZETA ( 3,X( l ) ,X( 2 > , ZB1 ) 

CALL ZETA { 3 , X < 3 I »X< 4 ) , Z B2 ) 

CALL ZETA ( 3 , X ( 5 ) , XI 6 > , Z A2 > 

TP(1,2>=ZB1(1)*R*X(1 )*X( 2) 

TP(1,4>=ZB2( 1>*R*X(3 >*XC 4) 

IF (LGFNITP(1,2)>X(2),K00E(3)»ZB1)) RETURN 
IF (LGFN(TP( 1,4) » X( 4 I , KODE (4)*ZB2) ) RETURN 
IF (LGFN(TP(1,3),X(&I, KODE ( 2 ) » Z A2 ) ) RETURN 
DSA2=CS(X<6> )-DLOG(DBLE(X<5) U-ZA2I4I 
DSBl=CSCX(2M-0L0G(DBLE(xm H-ZB114) 

DSB2=CS(X(4) I-DLOGI DBLE ( X( 3 } M-Z82C4) 

DHA2 = CH(X(6) ) + XI6)*( ZA2I l)~ZA2t 51 ) 

DHB1=CH( X ( 2 ) »+XI2>*l ZBim-ZBl(5M 
DHB2 = CH( X ( 4 ) I +X ( 4 ) *< ZB2 1 1 ) -Z B2 ( 5 )) 

Fill =DSB1-DSA1 
F(2I=DHBl-DHAl«-X{7)/2.0 
X13 = xm/XI3I 
SQX1 3=X13**2 

F(3J=TPCl,4)/R-TP(i, 2 )/R— X( l )*X( 7 ) *( 1.D-X13) 

F { 4) =DHB2-DHB1-X ( 7 ) * ( l.Q-SQX 131/2.0 
F (5»=DHB2-DHA2*X(7)*SQXl3/2.0 
F ( 6 I =DSB2-DSA2 

F(7)=ZA2I 1>*XI5)*X(6)-TP(1,3)/R 
IF ( MODI I ,5) .EQ.l) CALL RDFNS 
DO 7 J=1 , 7 
CONV ( J I =0 • 0 
DO 6 K=l»7 

6 CONVI J)=C0NV( J)+W( J,KI*F(KJ 

7 CONV < J l =CONV ( J ) /X ( J) 

TEST=0.0 

DO 8 J=i,7 

8 TEST=TEST+A6S (CONV ( J )) 

IF ( TEST • LT . EMAT ) GO TO 11 
DO 9 J=l, 7 

9 XI JI=XI JI*«1.0+CONV( Jl) 

10 CONTINUE 
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KODE ( 7 ) =1 

11 TP (4 ,2 ) =OHBl 

TP I 5 » 2 ) =OSBl 
TP(2,2)=X(1» 

TP(3,2»=X(2» 

CV=SNGLICP(TP|3,2) M-ZB1I6) 

GA=ZBl(3)+ZBl(2>**2/CV 

TP(7»2I=GA/ZB1(3) 

TP|8,2)=GA/ZB1<1) 

TP(6,2)=CV*TP(7,2» 

VM(3)=SQRT(X(7)*RI 

VM(1)=VM<3)/SQRT<TP<8,2»*TP< 1,2»/TP<2,2)) 
TP(3,3)=X(6) 

TP< 2,3>=X(5) 

TPI 4 , 3 ) =OHA2 
TP(5,3)=DSA2 

CV*SNGUCP(TP(3,3>n-ZA2(6> 

GA = ZA2I3M-ZA2(2)**2/CV 
TP(7t3)“GA/ZA2(3) 

TP(8,3)=GA/ZA2(1) 

TP(6,3)=CV*TP{7,3> 

TP(2,4)*X(3> 

TPI 3,4)=X<4) 

TPI 4,4)*DHB2 
TPI 5 ,4) =DSB2 

IF ITPI5,2).GT.TPC5,4)) KODE(B)=l 
CV*$NGLICPITPI3,4)n-ZB2l6> 
GA=ZB2|3)+ZB2I2)**2/CV 
TP(7,4)=GA/ZB2<3) 

TP ( 8 ,4I=GA/ZB2 ( 1 ) 

TP(6,4>=CV*TPI7»4) 

VMI6)=SQRTI2.D*R*(TP<4,3>-TPI4,4 ) ) ) 
VM|5)=VMI6)/SQRTITP(8»4)*TP| 1, 4I/TPI 2,41 > 
RETURN 
END 


DECK RGDFNS 

THE INVERSE OF THE DERIVATIVE MATRIX REQUIRED FOR DETERMINING THE 
VALUES OF X IS CALCULATED IN THIS SUBROUTINE. 

SUBROUTINE RDFNS 

COMMON /OUTNS/ E(47),ZA1(&),ZBK6),ZA2(6),ZB2(6), KODEI 8 ) 

COMMON /RONS/ W17,7> ,X(7 >,X13,SQX13 
DIMENSION XWI7I, YWt 7 I 
DOUBLE PRECISION CP 
DO 1 1=1,7 
DO 1 J=1 , 7 
1 W(I,J)=0.0 

CVB1=SNGL ( CP(X ( 2 ) > WB1I6) 

CVB2 = SNGL( CPI XI 4) ) I-ZB2I6) 

CVA2=$NGL ICPI X I 6 ) M-ZA2I6) 

W(1,1)=ZB1(2 )/X(l) 

W(l,2)*-CVBl/XI2) 

W|2,l)=< ZB1 I2)-ZB1(3 > J*X(2)/XU) 

WI2, 2)=-(CVBi+ZBll2) ) 

H I 2 ,7) =-0.5 

W(3,1)*X<2>*ZBH3>+XI7)*(1.0-2.3*X13> 
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W(3,2)=XI1>*ZB1<2> 

W(3, 3>=X(7)*SQX13-X( 4»*ZB2< 3) 

W(3,4)=-X(3)*ZB2I2I 

W(3»7)=XIU*{1.0-X13) 

W(4,t M-W(2,1)-XI7I*X13/X<3> 

W(4,2)=-W(2,2) 

M< 4,3>=<X<7 )*SQX13-X ( 4 > * ( ZB2 < 3 ) -ZB2 ( 2)) )/XC 3) 

W(4,4) =-( CVB2+ZB2 ( 2 ) > 

WI4,7)=(1.0-SQXl3l/2.0 

W(5»l)=W{2,l)+W(4,l) 

W ( 5 * 3 ) =W ( 4* 3 ) 

W ( 5 » 4) =W ( 4, 4 ) 

W(5»5)=X(6I*(ZA2(3)-ZA2(2))/X(5) 

W(5,M=CVA2 + ZA2<2> 

W(5,7)=W<4,7)-0.5 

W(6,3»=ZB2(2>/X<3> 

W(6,4>=-CVB2/X(4> 

W(6,5M-ZA2I2)/X<5> 

WI6,6) =CVA2/X(6) 

W(7,5)=-X(6)*ZA2(3) 

W(7»6)=-X(5)*ZA2(2) 

CALL MI NV (W,7,DET,XW,YW) 

MI NV IS THE SUBROUTINE USED TO INVERT THE DERIVATIVE MATRIX. THIS 
SUBROUTINE IS PART OF THE SCIENTIFIC SUBROUTINE PACKAGE SUPPLIED 

BY I.B.M. 

RETURN 

END 


DECK PGNS 

THE PERFECT GAS VALUES FOR THE MACH NUMBER, BOTH UPSTREAM AND 
DOWNSTREAM OF THE NORMAL SHOCK ARE CALCULATED IN THIS SUBROUTINE. 

SUBROUTINE PNS < R ATI 0, XM 1, XM2, DELR, KK ) 

COMMON /LDATA/ A(5»,G 

COMMON /GDATNS/ GL , G2 , G3 , G4, G5, G6, G7, G8, G9 , GIO,G 11 
COMMON /LIMITN/ EDAI , EMAT, EPNS 
DATA KG/O/ 

F< A) =S8*( A/( A+G5) ) **G4/ { G7* A-l . 0 ) **G3 
IF ( KG.NE.O ) GO TO I 
G1=IG-1. 01/2.0 
G2 = G 

G3=i.O/(G-l.O) 

G4=G*G3 
G5=2. 0*G3 
G6= { G+1 . 0 >*G3 
G7=2.0*G4 
G8=G6**G6 

G9=1.5*(G+l.0)**2/G 
G IO=F ( 4 . 0 1 
GI1=F( 5.0) 

KG= I 
1 R=RATI 0 
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KKK-0 
D R = U • D 

IF IR.LE.O.O.OR.R.GT.l.OJ GO TO 4 
IF (R.LT. 0.9999) GO TO 2 
X2 = 1.DHG9*( l.D-RH* *0.3 3333333 
SO TO 5 

2 Xl=4.3 
FI=G13-R 
X2=5.D 
F2=G11-R 

DO 3 I =1 » 2D 
0R=F2/R 

IF I ABS(OR) .LT.EPNS) GO TO 5 
DX=(X2-X1)/<F1/F2-1.0) 

X1 = X2 
F1 = F2 
X2= X2 + DX 

3 F2=FIX2)-R 
KKK = 2 

GO TO 5 

4 KKK=1 
X2-0.D 
Y2=0.D 
GO TO f> 

Y2= ( X2+G5 ) / ( G7*X2-L. 0 ) 
XMl=SQRT<X2) 

XM2=$0RT< Y2) 

DELR=DR 
KK=KKK 
RETURN 
ENO 


C DECK NSDATA 

BLOCK DATA 

COMMON /LIMITN/ E(3) 

DATA E/1.0E-6, 1.0E-5, l.OE-6/ 
END 


38 



APPENDIX D 


DESCRIPTION AND CARD LISTING OF THE COMPUTER 
ROUTINES THAT DESCRIBE THE VARIOUS GASES 

Any one of the sets of routines in this appendix is used with the routines in either 
appendix B or C. In addition, this set can also be used with the routines in reference 3 
whose deck names are RGASC1 and RDATA. This combination extends the applicability 
of reference 3 to include all the gases in this report. 

The deck names for these sets of routines consist of a one- or two-letter prefix 
identifying the gas and a root identifying the routine’s purpose. 


GENERAL PURPOSE OF ROUTINES 
Deck --ZETA 

In this subroutine, the compressibility -factor functions defined by equations (1) to (6) 
are calculated. With the exceptions of methane and natural gas, this subroutine can be 
used independently of the other routines. These points will be noted when the individual 
gases are discussed. Because of this independence, the use of this routine is described 
briefly here. This routine is referenced as follows: 

CALL ZETA(K,RHO,T,Z) 


where 

3 

RHO Density, kg/m 

T Temperature, K 

Z This is a six -element array. These elements are, in order, Zj to Zyj. 

K If K = 1, only Z(l) and Z(3) are calculated. If K = 2, only Z(2) and Z(4) are 

calculated. If K = 3, Z(l) to Z(6) are calculated. For the case of steam, 
Z(l) to Z(6) are calculated regardless of the value of K. 


Deck --TEMP 

The ideal-gas specific -heat functions defined by equations (7) to (9) are calculated 
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in this routine. Except for natural gas, whose ideal-gas specific heat is composition de- 
pendent, this routine can be used independently of the other routines. Because of this, 
the use of this routine is described briefly here. 

This is a double -precision function routine with three entry points - CP, CS, and 
CH. All three entry points use a single-precision argument, which is temperature. The 
function definitions are 


CP(T) - £ t (T) - Xi-ij 

(Dl) 

1 R 

s« 

CS(T) = 4 n (T) - id + In MU 
U R \RT/ 

(D2) 

H., 

CH(T) = £ m (T) = -i- - T 

(D3) 


From equations (8) and (9), it is seen that £ n and £ ni involve indefinite tempera- 
ture integrals. The constants of integration for these integrals are included in this rou- 
tine. The value of these constants is chosen such that the enthalpy and entropy equal zero 
at some reference thermodynamic state. Except for natural gas, steam, and low- 
temperature nitrogen, this reference state is that of the ideal gas at a temperature of 0 K 
and a pressure of 1x10^ pascals. The reference state for natural gas, steam, and low- 
temperature nitrogen will be given when these gases are discussed. 


Deck — LOG 


This is a logical function that tests whether or not the pressure and temperature lie 
within the range of the state equations. If the temperature range extends below critical, 
the vapor pressure relation is used to determine whether or not the fluid is a gas. Be- 
cause of the computational procedure, the lower pressure limit of the state equation for 
all gases has to be a nonzero positive number; this number is arbitrarily set at 0. 1 pas- 
cal. 


Deck --DATA 

This is a BLOCK DATA subprogram that supplies gas -dependent constants for the 
other routines. Among these constants are the molecular weight, the gas constant, and 
the perfect-gas specific-heat ratio. 
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Deck — TLG 


While the procedures in this report do not require this subroutine, the procedures 
in reference 3 do. The inclusion of this subroutine permits the application of the proce- 
dures in reference 3 to all the gases in this report. If necessary, this subroutine changes 
the temperature such that the fluid is a gas and lies within range of the state equation. 


ROUTINES APPLIED TO VARIOUS GASES 


Air 


The equations for the compressibility factor and the ideal-gas specific heat for air 
are the result of a curve fit of the data in reference 5. The temperature range of these 

5 

curve fits is from 150 to 1500 K. The pressure range extends to 100x10 pascals. These 
equations only apply to undissociated air. 

In order to check the accuracy of these fits, the values of Z and C p calculated by 
these routines were compared with those tabulated in reference 5. Over most of the pres- 
sure and temperature range, Z agrees to within 0. 1 percent and agrees to within 
0. 3 percent. The maximum disagreement is 0. 3 percent for Z and 3 percent for C^. 

The card listing of these routines follows. The prefix for the deck names is AI. 

C DECK AIZETA 

SUBROUTINE ZET A< KK, PP. TT ,Z ) 

DIMENSION Z(6I»A(10»3«2I»8(10t3t2) 

DATA A/0.0,9.742920BE-13,-4.2565612E-9,1.1837564E-7 t -1.19ll556E-6, 
14.2969317E-6,8.1108983E-6,-1.322802E-4,-1.402926E-4,1.194593E-3,5. 
2B580B55E-14,-6.044t752E-13,-1.2888247E-l3,4.6Q2D856E-l 1,-9. 90226 57 
3E-10,7.3369336E-9,-l . 9986894E-8,-5. 5490372E-8, 4. 8502964E-7 , 6. 55400 
416E-7.0.0, 2.4610094E-16,2.001808iE-14,-7.Q226616E-13,7.7383923t-12 
5,-3.0953497E-il,-2.6931749E-ll, 5. 104152E-10.-5.6194899E-10, I. 36719 
633E-10.0.0.-6.82004456E-12, 2 . 55393672E-8, -5.9187 82E-7. 4. 7646224E-6 
7, -1.28 90795 IE-5,-1. 622 17966E-5, l . 322802E-4, 0 . 0,1.1 94593E-3, -4. 6864 
8684E-l3 f 4.23092264E-12,7.7329482E-13#-Z. 3010428E - 10 . 3. 96090 628C-9 , 
9-2.35138008E-8,3.9973788E-8,5.5490372E-8,Q.0.6.5540Ql6E-7 f 0.0, -1.7 
A2270653E-15,-1.20108486E-13 f 3. 5 11 3308E- 12 . - 3. 095 35692E- 1 1 , 9.286049 
BlE-11.5. 3863 49 8E-11«— 5.1 04 152E-10»0.0»1. 367193 3E -10/ 

0ATA B/0.0»-7.79433664E-12,2.97959284E-8,-7.1025384E-7,5.955778E-6 
l»-1.71977268E-5,-2.43326949E-5, 2.645604E-4, 1 . 402 926E-4, 0. 0 , - 2. 63 61 
23848E-13.2.41767008E-12, 4. 5 1088645E-1 3 , - 1. 380625 68E-1 0 t 2. 475 56542E 
3-9,-1.5673B672E-8 f 2. 9980341E-8 » 5. 5490372E-8 , -2 . 4251482E-7 f 2*0.0, -6 
4.56269173E-16.-4.67088557E-14, 1 .40453232E-12 , - 1. 28973205E-1 1 , 4. 1 27 
513293E-ll,2.6931749E— il f -3.40276BE-lO, 1 . 8731633E-10 , 2*0 . 0 , 5 . 4560 35 
665E-11 » -1 • 7 877557 E-7 T 3. 5512692E-6 t -2. 38231 12E-5* 5.1 563180 4E -5.4. 86 
7653898£-5,-2.645604E-4,2*0.0,2.10891078E-l2,-1.69236906E-ll,-2.7D6 
853187E-12,6.9031284E-10,-9.9022657E-9,4.702160l6E-8,-5.9960682E-8 f 
9-5.5493372E-8,3*0.0,4.59388421E-15,2.80253134E-13,-7.0226616E-12,5 
A. 1589282E-11,-1.23813988£-10,-5.3863498E-11,3.4D2768E-10,2*0.0/ 
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K= KK 
P = PP 

$ = 1330. O/TT 

IFJK.EQ.2) GO TO 2 
BB = Ad, 1,1) 

CC= A( 1,2,1 ) 

DD= A| l ,3 1 1 ) 

00 1 N=2 , 10 

8B= BB*S+A(N,1,1) 

CC = CC*$+A(N,2, L) 

1 00= DD*S+A(N,3,1> 

2(11= l.O+(BB+(CC+DD*P)*P)*P 
Z ( 3 ) = 1.0+ ( 2.0*BB+< 3.0*CC+4.0*DD*P I *P )*P 
IF(K.EQ.l) RETURN 

2 BB= A(l, 1,2) 

CC= All, 2, 2) 

00= A(l, 3, 2) 

DO 3 M=2,10 

B8= BB*S+A( N, 1, 2 ) 

CC= CC*S+A(N,2,2> 

3 0D= DD*S+A( N, 3 , 2 ) 

2 ( 2 ) = 1. 0+( B8+( CC+D0*P ) *P) *P 
2 ( 4 ) = (B8+(CC/2.0+DD*P/3.0)*P)*P 
I F ( K.E3.2 ) RETURN 
BB= 611,1,1) 

BBB= B ( 1 , 1 , 2 ) 

CC= B(l,2,l) 

CCC= B (1, 2, 2 ) 

0D= Bd,3,l) 

DDD= Bd, 3, 2 ) 

00 4 N=2 ,10 

BB= BB*S+B ( N, 1 , 1 ) 

BBB= BBB*S+B ( N, 1, 2 ) 

CC= CC+S+B ( N, 2 , 1 ) 

CCC= CCC*S+B(N,2,2) 

DD= DD*S+B ( N , 3 , 1 ) 

4 DDD= DDD+S+B ( N, 3, 21 

Z ( 5 ) = ( BB + ( CC+DD*P > *P ) *P 
Z ( 6 ) = (BBB+(CCC+DDD*P)*P)*P 
RETURN 
END 


C DECK AITEMP 

DOUBLE PRECISION FUNCTION CP(T) 

DIMENSION A ( 8, 3, 2 ) 

DOUBLE PRECISION HI < 2 ) , SI ( 2 ) , S 

DATA A/0.0, 1.3466662 E-5, -3.1 114986E-4, 2. 41 2969 IE-3, -6. 332279E-3, 5. 
14472 78 3E-3, 1.4259098 E-3, 2. 4886014, 0.0,2. 24444367E-6, -6. 2229972 E-5, 
26.03242275E-4,-2.11075967E-3,2.72363915E-3,1.4259098E-3,2.48860l4, 
30.0,l.92383886E-6,-5. 185831 E-5, 4. 8259382E-4,-1.58306975E-3 , 1 . 81 5 75 
4943E-3,7 . 129549E-4, 2 .4886014»-l.752035E-8,6.39946R4E-7»-6.2254994E 

5- 6, 3. A 96617 4E- 5, -1.82616 13 E-3, 3. 2562223E-2 , -1 . 272201E-1 , 2. 62 5635 6 , 

6- 2.50290714E-9»1.06657807E-7»-l.24508188E-6»8.7415435E-6, -6.087204 
733E-4, 1.6281 11 15E-2, -1.2 72201E-1, 2. 6256356, -2. 19004375E-9, 9. 142097 
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87lE-8,-1.03>56823E-6,6.9932348E-6,-4.56540325£-4, 1 . 08540743E-2 . - 6 . 
93610056-2,2.6256356/ 

DATA HI, SI/-. 8930568474993379, - 8. 3696065 5648124*, ,21. 32793295150501 
1,21.43350926781483/ 

L = 1 

1 S=T/1. 0D2 
K = 1 

l Ft S.GT.5.12397093D0 )K=2 
CP=A(1,L,K> 

00 2 N=2,7 

2 CP=CP*S+AIN,L,K) 

GO TO (3,4,5),L 

3 CP=CP*S+A(8,l,K) 

RETURN 

ENTRY CSIT) 

L a 2 

GO TO l 

4 CP=CP*S+A{B,2,K) *DLOG f S I *S I f K ) 

RETURN 

ENTRY CHIT) 

L = 3 

GO TO 1 

5 CP=T*(CP*S+AI8,3,K)I+HI(K> 

RETURN 

ENO 


C DECK AILOG 

LOGICAL FUNCTION LGFN <P,T,J,Z) 

DIMENSION Z ( 6 ) 

J = 1 

LGFN=. TRUE. 

IF (P.LT.0.1.0R.P.GT.1.1E7.0R.T.LT. 149. O.OR.T.GT .1510.0. OR. Z(1).LE 
l.O.0.OR.Z(2).LE.0.0.OR.Z(3).LE.0.O) RETURN 
J=0 

LGFN=. FALSE. 

RETURN 

END 


C DECK AIDAT A I 

BLOCK DATA 
COMMON /LDATA/ RI6) 

DATA R/l. 0,287. 0588, 28. 9641, 1.6333146-3,5.6,1.4/ 
END 
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C DECK AITLG 


SUBROUTINE TLOGIC(P,T) 

IF(T.LT.149.0)T=149.0 

RETURN 

END 


Nitrogen 


There are two sets of nitrogen routines: one covers the cryogenic temperature 
range, and the other covers the high-temperature range. 

Nitrogen, range I . - The equations for the compressibility factor and the ideal-gas 
specific heat for nitrogen (range I) are from reference 6. This reference claims a tem- 
perature range from 64 to 350 K and a pressure range to 300xl0 5 pascals. This report 
assumes a temperature range from 60 to 400 K and a pressure range to 300xl0 5 pascals. 
These are the same ranges that were assumed in reference 7. 

The values of Z and C p calculated by these routines (using the state equations of 
ref. 6) were compared with those tabulated in reference 5. Over most of the pressure 
and temperature range where comparisons are possible, Z agrees to within 0. 2 percent 
and Cp agrees to within 1 percent. The maximum disagreement was 0.4 percent for Z 

and 2 percent for C . 

P 

The reference state for the enthalpy and entropy is that of the saturated liquid at the 
triple point. 

The card listing of these routines follows. The prefix for the deck names is N. 

C DECK NZETA 

SUBROUTINE ZETAI KK,PP,TT»Z ) 

DIMENSION Z ( 6 ) 

DAT A A,Bl,B2,B3»B4,B5»B6,B7,B8,89,B10,Bll,Bi2,B13,B14,B15/-7.135E- 
16,1.2D34917E-3,-2.510789E-l,-4.9681584El,3.7073373E2»1.496473E6, 
22.1027719E-6,-2.4516046E-4, 2.3102822E-9, 4.9866482, 1.6771 286E3, 
3-l.656225E5,-6.5374809E-5, 2. 42091 088-2,-1.126389, 1. 1829604E-12/ 

C0( X, Y,Z)=X*S3+Y*S4*Z4S5 
K=KK 
P = PP 
T=TT 

Sl=1.0/T 
S2=S1*S1 
S 3=S1*S2 
S4=S1*$3 
S5=S1*S4 
P2=P*P 

EXPO=EXP| A*P2) 

AA=1.0/(2.D*A) 

BB= 1 • D/A 

C1=C0( B9, B1 0, 81 1 ) 

C2=C0(B12,B13,B14) 
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IF (K.EQ.2) GO TO 1 

BA1=B1+B2*S1+B3*$2+B4*S3+B5*S5 

BA2=B6+B7*S1 

ZA=1.0+(BAI+(BA2+(BS+815*S1*P2)*P>*P)*P 
ZB=EXP0*P2*{Cl+C2*P2 > 

Z<1)= ZA+ZB 

ZA=1.0+(2.0*BAl+(3.0*BA2+< 4.0+B8+6. 0*B15*S 1*P2 ) *P 1 *P 1 *P 
ZB=P2*EXPO*< 3.0+C1 + I 2.0*A*C1+5.0*C2+2.0*C2*P2*A)*P2) 

Z ( 3 1 = ZA+ZB 
IF (K.EQ.l) RETURN 

1 BT=B1-B3*S2-2.0*B4*S3-4.Q*B5*$5 

C1P = -CO(3.0*B9,4.Q*B10,5.0*BU) 

C2P=-C0(3.0*B12,4.0*B13,5.0*B14) 

CC1=C1+C1P 

CC2=C2+C2P 

ZA=1.0+(BT+(B6+BB*P) *P 1 *P 
ZB=P2*EXPO*(CCl+CC2*P2> 

Z 1 2 ) - ZA+ZB 

ZA= (BT+( 86/2. 0+B8*P/ 3.0) *P )*P 
ZB=AA*IEXP0*(CC1+{P2-BB)*CC2)+BB*CC2-CC1 ) 

Z ( 4 1 * ZA+ZB 
IF (K.EQ.2) RETURN 
C1PP=C0(9.D*B9, 16.0+610,25.0*6X1) 

C2PP=C0I9. 0*812,16. 0*B13, 25. Q*B14) 

ZA=-I(B2*Sl+2.0*B3*S2+3.0*B4*S3+5.0*B5*S5)+( B7*S 1/2. O+Bl 5*S1*P**3/ 
15.0 1 *P 1 *P 

ZB- A A* ( EXPO* (C1P+(P2-BB)*C2P)+BB*C2P-C1P) 

Z 1 5 > - ZA+ZB 

CC1P=C1P+C1PP 

CC2P=C2P+C2PP 

ZA= ( 2. 3*B3*S2+6.0*B4*S3+20.0*B5*S5 1 *P 
ZB=AA*(EXP0*{CC1P+(P2-BB)*CC2P)+CC2P*BB-CCIP) 

Z ( 6 I = ZA+ZB 

RETURN 

END 


C DECK NTEMP 

DOUBLE PRECISION FUNCTION CP(T) 

DOUBLE PRECISION A<5,3>, S,HI,SI 

DATA A, SI, HI /2. 50114600, -9. 720561D-5, 1. 0360560-6, -4.437258D-9, 
16, 82 55 96D-1 2 , 2.501146 DO » -9 .720581 D- 5,5.18028 D- 7,-1.4790860-9, 
21.7063?90-12,2.50114600,-4.8602905D-5,3.45352D-7,-l. 1093145D-9, 
31.36511920-12,7.71240-1,5.083102/ 

K= 1 

1 S=T 

CP = A ( 5 , K) 

DO 2 1=1,3 

M= 5-1 

2 CP= CP*S+A1 M»K ) 

GO TO ( 3,4, 5 1 ,K 

3 CP= CP*S+A< 1,1) 

RETURN 

ENTRY CS ( T 1 
K=2 
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GO TO 1 

4 CP= CP*S+A(l, 2 )*DLOG(S)+SI 
RETURN 

ENTRY CH( T ) 

K = 3 

GO TO 1 

5 CP= (CP*S + A(l,3) I *S *H I 
RETURN 

END 


C DECK NLOG 

LOGICAL FUNCTION LGFN (P»T,J,Z> 

COMMON /LDATA/ XKV,R , XMW , RC, D2, G 
DIMENSION Z ( 6 ) 

S = T 
J = 1 

LGFN=. TRUE. 

IF < P.LT.0.l.OR.P.GT.3.1E7.0R.S.GT.410.0.0R.S.LT.59.0.0R.Z< D.LE.O 
l.O.OR.Z<2).LE.0.0.OR.Z(3I.LE.0.0l RETURN 
IF < S.GT. 126.261 GO TO 1 

PS*10.D**(-305.07339/S+5. 5335216+1 . 16441 101+ I -3. 1 3B92Q5E-3+ < 2. 9857 
1103E-5+l~1.4238458E-7+2.7375282E-l0*S»*S)*SI*S)*S) 

IF ( P. GT. XK V*PS I RETURN 
1 J=0 

LGFN*. FALSE. 

RETURN 

ENO 


C DECK NDATA 

BLOCK DATA 
COMMON /LDATA/ RI 6 ) 

DATA R/ 1.0, 296.8008, 28.0134, 1. 579703E-3, 5. 6 , 1.4/ 
END 


C DECK NTLG 

SUBROUTINE TLOGIClP,T) 

DIMENSION A ( T > 

DATA A/5.8435l56E-6 t -9. 652 1356E-5, 3. 2061 122E-3 8. 1845679E-2 , 
11.2169369,-5.2033603,50.343493/ 

IFIT.GT.126.26) RETURN 
V= ALOG(P) 

S* A(ll 
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DO 1 N=2 , 7 

1 S= S*V*A(N) 

IF(T.L T .S) T-S 

IF( T.LT.59.0) T*59.0 

RETURN 

END 


Nitrogen, range n . - The equations for the compressibility factor and the ideal-gas 
specific heat for nitrogen (range II) are the result of a curve fit of the data in reference 5. 
The temperature ranges of these fits is from 170 to 3000 K. The pressure range extends 
to HOxlO 5 pascals. These equations only apply to undissociated nitrogen. 

In order to check the accuracy of these curve fits, the values of Z and Cp calcu- 
lated by these routines were compared with the values tabulated in reference 5. Over 
most of the range, Z agrees to within 0. 1 percent, and Cp agrees to within 0. 5 percent. 
The maximum disagreement is 0. 7 percent for Z and 2 percent for Cp. 

The card listing of these routines follows. The prefix for the deck names is HN. 

C DECK HNZETA 

SUBROUTINE ZETA <KK,PP,TT,Z) 

DIMENSION BB(11,4), CC(il,4), 00(11*4), BCD(11,4,3>, Z16>, A(4,3) 
EQUIVALENCE ( BB< 1 ) ,BCD( 1 > ) , ( CC< 1 ) , BCDI 45 ) ) , ( DD(1) ,BCD( 89) > 

DATA BB/D.,1. 338343E-9 ,-6. 4553703E-9, -1 . 38D7339E-7, 1. 0851228E-6, -1 
1. 41085 47E-7.-1. 85851 43E-5, 3 . 015 3933E-5, -4. 9439856E-4, 1 . 6443824E-3 , 
2-1.0717843E-4,0.,-9.3684DlE-9,3.B732225E-8,6.9036695E-7,-4.3404912 
3E-6,4.2325641E-7»3.7170286E-5,-3.0153933E-5,0.,1.6440324E-3,-2.143 
45686E-4,Q»*-1, 07Q6744E-8 *4, 5187596E-8»8. 2844034E-7, -5. 425614E-6*5. 
564341BBE-7,5.5755429E-5,-6.0307866E-5,4.9439856E-4,0.,-1.0717843E- 
64,0, ,7.4947208E-8,-2.7112557E-7,-4.1422017E-6,2.17t>2456E-5,-1.6930 
7256E-5,-1.1151086E-4,6.0307866E-5, 2*0. , -2, 1 435686E-4/ 

DATA CC/0.,-1.4161135E-li,-1.8207963E-10,-9.l721547E-lQ,5.34tl861E 
1-8, -2. 000952E-7,-8.9753116E-7,6.4975889E-6,-l. 1598976E-5, 1.065 3882 
2E-5, -2. 7372964E-6, 0., 0. 49668 IE-1 1, 9.1 039815E-10, 3. 668861 9E-9, -1.60 
323558E-7,4.001904E-7,8.97531l6E-7, 0.,-1.1598976E-5,2.1307764F-5,-8 
4. 2118B92E -6,0. ,4.956 3972 E- 11, 5. 4623889E-10»2.2930387E-9, -1.0682372 
5E-7 , 3. 001428E-7, 8. 9753 11 6E-7 , -3. 2487944E-6, 0. , 5. 326941E-6, -2. 73729 
664E-6,0.,-2.9738384E-10,-2.73ll944E-9,-9.1721547E-9,3.2047117E-7,- 
76.002856E-7,-8.9753116E-7, 2*0. , 1. 0653882E-5, -8.21 18892E-6/ 

DATA DD/8.9688983E-14, 3. 7837236E-13, -6. 3374971E- 12,4.2678663E-ll,- 
19.4492282E-lO,3.5667223E-9,2.022306BE-8,-1.3522968E-7,2.6433l66E-7 
2,-2.O302481E-7,4.9983348E-8,-5.381339E-13,-1.8918618E-12,2.5349988 
3E-11,-1.2803599E-10, l. 889845 6E-9* -3. 566722 3E-9,0. ,-1.352 2968E-7, 5. 
42866332E-7,-6.1147443E-7,1.9993339E-7,-2.0927429E-13,-7.567447E-13 
5, 1.0562495E-ll,-5.6904884E-ll,9.4492282E-10,-2.3778149E-9,-6.74i02 
627E-9, 0., 8. 81 10553E-8, -1. 358832 IE-7, 4. 99833486-8,1. 2556458E-1 2, 3. 7 
78 372 36 E-l 2 ,-4.224998 IE-1 1, 1 . 7071465E-10 , -1. 8898456E-9 , 2. 3778149E -9 
8,2*0., 1.76221 11 E-7,-4.0764962E-7, 1. 9993339E-7/ 

S-1000.0/TT 

P=PP/S 

K=KK 

DO 3 J=1 ,4 
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IF (K.EQ.2.AND.J.NE.2) GO TO 3 
DO 2 1=1,3 
A I J, I ) =BCDC 1 , J, I ) 

00 1 H=2 til 

1 a ( j, n =A( j, n*s*BCDi m, j, n 

2 CONTINUE 

IF IK.EQ.l.AND.J.EQ.l) GO TO A 

3 CONTINUE 

Z(2)=1.0+(AI2,1)+(A(2,2)+A(2,3)*P)*P)*P 
Z (4)=I AI2»i )♦! A(2,2)/2.0+A(2»3)*P/3.0)*P)*P 
IF IK.EQ.2) RETURN 
ZI5)=(A(3,1)+(AI3,2)+A(3,3)*P)*P)*P 
Z(6)=<A{4,1]+<A(4,2)+A{4,3)*P)*P)*P 

4 Zm=1.0+(A(l,l)+(Af 1,2)+AU,3)*P)*P)*P 

ZI3)=1.0'M2.0*A(1,1)M3.Q*A(1»2)+4.Q*AI1,3)*P)*P)*P 
RETURN 

ENO 


C DECK HNTEHP 


1 


2 

3 


4 


5 


DOUBLE PRECISION FUNCTION CP IT) 

DIMENSION A ( 8 , 3, 2 ) 

DOUBLE PRECISION CS» CH, S ,H 1 1 2 ) , SII 2 ) 

DATA A / 19. 2 04975 » -43. 37888, 31.58 32 31,-6.7699 321»-5.2724127E-3, 6. 42 
1635 6 7E-2 f 1 .43 1 144E-2 t 2. 4986544, 2. 7435679, -7. 2298135, 6. 3166462, -1.6 
292483, -1. 75 74709 E- 3, 3. 21 31 78 4E -2, 1 • 431 144E-2 ,2 .4936544,2.40062 1 9*- 
36. 1969 B3, 5. 2638718,-1. 3539864, - 1. 3181032E-3, 2. 1421189E-2 , 7. 15572 E- 
43,2. 4986544, l. 227 4093E-3, -6.902 4646E-4, -9. 369791 4E -2, .5882261, -1.4 
5746878, 1.4646469,.29 16646, 2. 1558557, 1.7534419E-4,-1.1504108E-4,-l. 
68739583E-2,. 14705652, -.4915626,. 73232345,. 2916646, 2. 1558557, 1.5342 

I }i!^;; 9 ; 8 ™c^: 5, " l ’ 56163l9E_2t - 117 ^ 522 »-- 36867i95 »^ 882 i563, 

“• l Zm 1558557/ 

DATA HI, SI/-. 923621 54000 1477 1,1 00. 93794267 75437, 2 6. 17073912726212, 
125.67752552153525/ 

K = 1 
J = 1 


S=T/1. 3D3 

IF IS. GT.. 8119820835) J=2 
CP=A ( 1 ,K, J ) 

DO 2 N=2,7 
CP=CP*S+AIN,K, J) 

GO TO I 3, 4, 5 ) , K 
CP=CP*S+A 1 8 , 1 , J ) 

RETURN 
ENTRY CS(T) 

K=2 


GO TO 1 

CP = CP*S+A(8,2, J)*DLOGIS)+SH J) 

RETURN 

ENTRY CHIT) 

K=3 


GO TO 1 

CP=T*(CP*S*AI8,3,J))*HII J) 

RETURN 

END 
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C DECK HNLOG 


LOGICAL FUNCTION LGFN (P,T,J,Z) 

DIMENSION Z ( 6 ) 

J = 1 

IF IP. LT. 0. 1 .OR.P.GT . 1.1E7.0R.T.LT.1 69.0. OR. T.GT . 3010.0. OR. Z(1)*LE 
1.0.0.0R.ZI2) .LE.O.O.OR.Z(3).LE.O.OI RETURN 
J =0 

LGFN=. FALSE. 

RETURN 

END 


C DECK HNDATA 

BLOCK DATA 

COMMON /LDATA/ R C 6» I , , , ,, 

DATA R/l.Ot 296. 8008. 2 8. 013 A, 1. 579703E-3, 5. 6 f 

END 


C DECK HNTLG 

SUBROUTINE TLOGIC (P,T) 
IF IT. LT. 169.0) T=169.0 
RETURN 
END 


Oxygen 

There are two sets of routines for oxygen: one covers the cryogenic temperature 
range, and the other covers the high-temperature range. 

Oxygen, range I . - The equation for the compressibility factor for oxygen (range I) 
is from reference 8, and the equation for the ideal-gas specific heat is the result of a 
curve fit of the data in reference 5. Reference 8 claims a temperature range from 65 to 
300 K and a pressure range to 300xl0 5 pascals. This report assumes a temperature 
range from 60 to 400 K and a pressure range to 300x10^ pascals. 

The values of Z and C p calculated by the routines in this report (using the state 
equation of ref. 8) were compared with those tabulated in reference 5. Over most of the 
pressure and temperature range where comparisons can be made, Z agrees to within 
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0.2 percent and C p agrees to within 1 percent. The maximum disagreement is 3 per- 
cent for Z and 3 percent for C p ; the maximum disagreement in Z occurs at 170 K 
and 70 bars. 

The card listing of these routines follows. The prefix for the deck names is O- 

C DECK OZETA 


SUBROUTINE ZET A < KK, PP t TT t Z ) 

DIMENSION Z<6> , A < 3, 3 , 4 1 , FI ( 3 > , F2 ( 3 ) 

DATA A1,A2, A3.A32, A35 , A36, A4, A42, A45, A46, A5, A52, A55 , A56, A6, A62, A65 
l»A7,A8,A85,A9,A10,Al02,A105,A106,All,A12,A125,Ai3,A135»A14,A145,A2 
24,A244,A25, A26, A27, A29, A28 , A30, A31 , ALPHA, TC/1. 29024534E- 3 , -.5012 53 
3921 ,-2814.00595,5628.01189, 8442.01784,-16884.0357, 7314650.15,-2925 

48600. 6, -365 73250. 7, 14629 3003., -1.10552497E 10, 6. 6 3314981E 10, 7.73867 

5478E10,-4.64320487EU,-6.78596713E-10,-1.3571934 3E-9,-3.392983 57E- 

610,9.48465081E-7,7.22544527E-5,-3.61272264E-5,-3.225972l8E-2,-. 42 7 

78 204 94,. 8 556 40987,. 64173 07 4, -1.283 46 148, 3. 8 02 1 54 84E-10, 7. 0848317 9E 

3616106E-8,1.41493762E-10,-3.53734405E-ll,2.83963097E-8,-1.4 

91981549E-8,-5.46946646E-6,-1.09389329E-5,2.37789478E-5,-3.0667678E 

*9,6.02168773E-8,3.01534684E-2,232.8 32966,154.77/ 

DATA A/2.06703766,3591.79212,-416047.043,1.03087315E-5,-3.113551E- 

13,1.22849158,6.40l42244E-li,-1.27186613E-8,1.66719727E-6,-4.134075 

231,-13775.3764,1664l88.17,-2.06174629E-5,9.34065301E-3,-4.91396632 

3 * -1 • 28 32B449E- ID, 3. 81 559 83 9E-8, -6. 668789 09 E-6 ,-6.20111297 ,-143 67.1 

22 3.09261944E-5, 1 . 2454204E-2 , - 6.14245 79, -1.9 2 042673E 

f~i t> , t !=?!I 46452E_8,_8#33598636E " 6,12 * 4022Z59 » 4 ^101 .5054,-8320940. 86 

7^ 6 T l 2 5 >!!«5lh‘!*, 7362612E " 2,24 * 56903l6,3 * 84O85346E-lO »- l « 526 23935 

' c-7 t 3* 334394558-5/ 

AF ( P )= ( Bl+I B2+( B3+84*P)*P) *P)*P 

BF1 ( A, B ,C ) - ( A+< B+C*P2 »*P2) *P2*EXPB 

BF2 ( A, B,C )=A*S3+B*S4+C*S5 

BF3( A, B,C,D)=< (A+(B+C*P2>*P2-(B+D*P2-D/A24)/A24)*EXPB-A* (B-D/A24)/ 
1A24 ) /A244 
K=KK 

p 3 pp 

T = TT 


S=1 .0/ T 

S2=$*S 

S3=S2*S 

S4=S3*S 

S5=S4*S 

S7=$5*S2 

P2=P*P 

APB=A24*P2 

EXPB=EXP(APB) 

P28=P**A28 
PALPHA=P28- ALPHA 
APC= 42 S*PAL PHA**2 
EXPPC=EXP«APC) 

TTC=T-TC 
ATC=A27*TTC+*2 
EXPTC = £XP< ATC ) 

ZC1 = A2 5*S*P28*PALPHA*EXPPC*EXPTC 

I F ( K.EQ.2 ) GO TO 2 

81=A1 t A2*S+A3*S3+A4*S5+A5*S7 

B2 = A6*T«-A7 + A8*S + A9*S2+A10*S3 

B3=A11 f A1 2*S 

B4=A13*S*A14*S2 
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ZA=1.D+AF(P) 

00 1 N = l» 3 

1 F1(N)=BF2(A(1,N, 1),A(2,N,1)#A<3,N,1) > 

ZB = BF1(FH 1 )tFl(2)»FlC3) ) 

Z ( 1 I =ZA+ZB+ZC1 
IF(K.EQ.O) RETURN 
Bl-2. 3*B1 
82=3. 0*B2 
B3=4.D*B3 
84=5.0*84 
ZA=AF( P1+1.0 

ZB= ( l.0+A244*P2l *ZB+BF1( 2.0+F1 ( L) » 4.0*F1 ( 2 ) * 6.0*F1< 3 ) ) 
ZC=< l.0+A28*(1.0+P28*< 1. 0+2. 0*APC ) /P ALPHA ) )*ZC1 
Z(3) =ZA+ZB+ZC 
I F ( K.EQ. 1 (RETURN 

2 81=A1+A32*S3+A42*S5+A52*S7 
B2=A62*T+A7-A9*S2+A102*S3 
B3=A11 

B4=-A14*S2 
ZA=1 . D+AF (PI 
DO 3 N=l»3 

3 Fl(N>=BF2(A(l f N f 2) ,A< 2,N*2 >, A( 3,N,2) ) 
ZB=BF1(F1(1)»F1(2)»F1(3) > 

ZC=A29*T*TTC*ZCi 

Z < 2 ) =Z A+ZB+ZC 

B2=B2/2.0 

63=83/3.0 

B4=B4/4.0 

ZA=AF(P) 

ZB=BF3(F1 (1 ),FlC2)tFl(3) ,2.0*F1( 3) ) 
ZC=A3L*TTC*EXPTC*(EXPPC-A30) 

Z { 4 ) =Z A+ZB+ZC 

IF(K.EQ.2»RETURN 

B1=A36*S3+A46*S5+A56*S7 

82= A6*T+A9*S2+A1Q6*S 3 

B3=0.0 

B4=-Al45*S2 

Z A=AF( P) 

DO 4 N*l,3 

Fl(N)=BF2CA(l,N,3)*A(2,N,3J*A(3*N,3n 

4 F2<Nl=SF2(A(l,N,4»,A(2,N,4),A<3,N f 4) > 
ZB=BF3(F2(1»,F2(2»,F2(3),2.0*F2(3) ) 

ZC = T*Z;*( 1. £>+2.0*ATC ) /TTC 

Z (6 »=ZA+ZB+ZC 

B1=-A2*S+A35*S3+A45*S5+A55*S7 

B2=A65*T+A85*S-A9*S2+A105*S3 

B3=A125*S 

B4=A135*S+A145*S2 

ZA=AF( P) 

ZB=BF3(F1(1 )«Fl(2)»FL(3) »2.0*Fi ( 3 1 ) 
ZC=A31*(TTC-S/A29)*EXPTC*< EXPPC-A30 ) 

Z (5 ) = Z A+ZB+ZC 

RETURN 

END 



C DECK OTEMP 


DOUBLE PRECISION FUNCTION CP(T) 

D0U8LE PRECISION A( 9 , 3 > , HI , S I , S 

DATA A»HI ,SI/3. 05760850-5, -3. 62300430-4, 1.3986432D-3, -1.82346020-3 
1, 2. 57413350-3, -1. 3617530-2, 3. 25342050-2, -3. 2270317D-2, 2. 512827300, 
23. 8220 10620-6, -5. 175 72043D-5, 2. 331Q72D-4, -3. 64692040-4, 6. 43533375D 
3-4, -4. 53917 6 67 0-3, 1.62671025 0-2,-3.22703170- 2, 2. 512 82 7 30 0,3. 3973 ^2 
4780-6,-4.528755380-5, 1.998061710-4,-3.039100330-4,5.1482670-4,-3.4 
50438250-3, 1.0844735D-2, -1.61 35 15850-2, 2. 5128273D0, -1.515548600, 22. 
620265612400/ 

K=1 

1 S=T/ 100.000 
CP=A(1,K) 

DO 2 N=2, 8 

2 CP=CP*S+A(N»K) 

GO TO 13,4, 51, K 

3 CP=CP*S+A(9,i) 

RETURN 

ENTRY CS( T) 

K=2 

GO TO 1 

4 CP-CP+S+A(9,2)*DL0G( S )+Sl 
RETURN 

ENTRY CHIT) 

K=3 

GO TO 1 

5 CP=T*( CP*S+A (9, 3 ) I+HI 
RETURN 

END 


C DECK OLOG 

LOGICAL FUNCTION LGFN <P,T,J,Z) 

COMMON /LDAT A/ XKV.R ,XMW ,RC, 02, G 
DIMENSION 7(6), A(9) 

DATA A/. 5 115947 ,-1.469853, -1.8221345, 11.180811,-13.613656,3.313939 
13,3.3609539,19.85226,-8.8689461/ 

S=T/100.0 

J*1 

LGFN=. TRUE. 

IF (P, LT.O.l.OR.P.GT .3. 1E7.0R.S.LT. 0.59. OR. S.GT. 4. l.OR. Z ( 1 > . LE. 0.0 
1. OR. Z< 2).LE.O. O.OR.Z I3).LE.0.0) RETURN 
IF < S.GT. 1.5477) GO TO 2 
X=A ( l ) 

DO 1 N=2, 8 

1 X=X*S*A(N> 

X=EXP( X+A (9 ) /S ) 

IF ( P. ST. XKV*X ) RETURN 

2 J = 0 

LGFN=. FALSE. 

RETURN 

END 
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C DECK ODATA 


BLOCK DATA 
COMMON /LDATA/ R(6I 

DATA Rfl. 0,259. 8347, 31.9988, 1. 804444E-3* 5. 6, l . 4/ 
END 


C DECK OTLG 


SUBROUTINE TLOGICIP.T) 

DIMENSION AI8> 

DATA A/l. 5415673E-7* -5.58328E-6* 2. 7328988E-4»-3. 4963316E-3,- 
l!)68E-4,5.487465E-lt -1.4597974, 49. 35502/ 

I F(T.GT.154.77)RETURN 
V-ALOGI P) 

S = A ( 1 ) 

DO I N=2 » 8 


1.5388 


1 S=S*V*A(N> 

IF(T.LT.S) T=S 
IFIT.LT.59.0)T=59.D 


RETURN 

END 


Oxygen, range II . - The equations for the compressibility factor and the ideal-gas 
specific heat for oxygen (range II) are the result of a curve fit of the data in reference 5. 
The temperature range of these fits is from 180 to 3000 K. The pressure range extends 
to 100x10^ pascals. These equations only apply to undissociated oxygen. 

In order to check the accuracy of these curve fits, the values of Z and C p calcu- 
lated by these routines were compared with those tabulated in reference 5 . Over most 
of the temperature and pressure range, Z agrees to within 0. 1 percent and C p agrees 
to within 0. 5 percent. The maximum disagreement is 0. 3 percent for Z and 16 percent 

for C . This 16 percent disagreement in C occurred at a temperature of 180 K and a 
P P 

pressure of 40 bar. 

The card listing of these routines follows. The prefix for the deck names is HO. 

C DECK HOZETA 


SUBROUTINE ZETA IKK,PP,TT,Z» 

DIMENSION BBU1,4), CC(11,4), DD(lt,4J, BCD(11,4,3», Z<6), A(4,3I 
EQUIVALENCE (BB( 1 ) * BCD( 11)* (CCI 1 ) , BCD( 45)) , (D3(1),BCD(89»> 

DATA BB/-6.0201768E-10 , 3.97764E-9, 3. 168l402E-8»-2. 5592444E-7,-8. 15 
113199E-8,7.9876699E-7,2.4938354E-5,-1.6164074E-4 f -1.606299E-4,l. D6 
2l436E-3»-2.8220151E-5, 4.816141 4£ -9, -2.784 348 E-8,-1. 90088 41E-7 » 1.27 
396222E-6,3.260528E-7,-2.396301£-6,-4.9876708E-5, 1.6164074E-4,D. , l. 
4061436E-3 » -5.6440302 E-5,5.418l591E-9,-3*182il2E-8»-2.217698lE-7»l. 
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553554566-6, 4.G7566E-7 , -3 . 195068E-6, -7. 4815062E-5 , 3. 232 8148E-4 , 1 . 60 
662996-4,0.,-2.8220151E-5,-4.3345273E-8,2.2274784E-7,l.3306189E-6,- 
77. 6777332E-6.-1.630264E-6, 9. 58520396-6,1. 4963012E-4, -3. 2328148E-4, 
82*0. » -5. 6440 302 E -5/ 

DATA CC/2*0.,3.00822 39E-ll,-3.43819UE-9,3.3326542E-8,-9.0112977E- 
18,3.5646415E-7,-2.4130715E-6»6« 4195294E-6, -4. 022482 IE-6, 9. 91 2 1 37E- 
27,2*0. , -1. 50411 19E-10, 1. 3752 764E-8, -9. 9979626E-B, 1.80225956-7,-3.5 
3646415 E-7 , 0., 6. 41952 94E-6, -8. 0449642E-6, 2. 973641 IE-6, 2*0., -9. 02467 
41 7 E- 11 ,8.595477BE-9, -6. 665 3084E-8 , 1 . 351 6947E-7 , -3. 564641 5E-7 , 1 . 2 06 
55358E-5,0.,-2.0U241E-6,9.912137E-7,2*0.,4.512 335BE-10,-3.4381911E 
6-8,1.9995925E-7,-2.7033893E-7,3.5646415E-7,2*0.,-4.0224821E-6,2.97 
736411E-6/ 

DATA D0/3*0., 1.6667272E- 10, -1.9308409E-9, 5. 48680 828-9,3. 0655349E -9 
1,-7. 5825703E-9,-5*6411225E-8, 8.3326385E-8, -2. 6382638F-8 ,3*0. ,-5.00 
201816E-10,3.8616818E-9,-5.4868082E-9,0.,-7.5825703E-9,-1.1282245E- 
37,2.4997916E-7,-1.0553055E-7,3*0.,-2.2223029E-10, 1. 9303409E-9, -3.6 
4578721 E-9,-1.021845E- 9, 0., -1.8 80374 2E -8 ,5.5550923E-8,-2.6382638E-8 
5,3*0.,5.6669088E-10,-3.8616818E-9, 3.6578721E-9, 2*0. ,-3.7607483E-8, 
61. 6665277E-7, -1.05530556-7/ 

S=1003.0/TT 

P=PP/S 

K=KK 

DO 3 J = l,4 

IF (K.EQ.2.AND.J.NE.2) GO TO 3 
DO 2 1=1,3 
A(J,I)«BCD(l,J,n 
DO 1 M=2 ,11 

1 Atj,n=A(j,n*s+BCO(M,j, n 

2 CONTIMOE 

IF ( K. EO. 1. AND. J.EQ. 1 ) GO TO 4 

3 CONTI NJE 

Z (2) =1.0 + 1 A < 2, 1 ) + < A( 2,2>+A( 2,3I*PJ*P)*P 
Z<4)=( A(2,l > + < A (2, 2) /2.0+A<2,3»*P/3.Q)*P)*P 
IF ( K. EQ* 2 1 RETURN 

Z{5)=(A(3,l)+(A(3»2)+A(3»3)*PI*Pl*P 
Z<6)=< A< 4,11+1 A( 4,2) +A< 4,3 )*P)*P)*P 

4 Z<l)=l.0+IA(1,1)+<A( l,2)+A(l,3)*P)*P)*P 
Z(3)=1.0+(2.0*A(1»1) +( 3«0*A( 1, 2 1+4.0* A 1 l, 3) *P )*P )*P 
RETURN 

END 


C DECK HOTEMP 

DOUBLE PRECISION FUNCTION CP <T) 

DOUBLE PRECISION HI(2),SI(2),S 
DIMENSION A ( 9, 3, 2 ) 

DATA A/-223. 15231,414.9894,-195.44586,-64.414824,67.444744,-5.4398 
1936,-3.9799316,0.8979885,2.4467825,-27.894039,59.2842,-32.57431,-1 
22.882965, 16 . 861186,-1.8132979, -1. 9899658 , . 8979885,2.4467825,-24. 79 
34701,51.873675,-27.920837,-10.735804,13.488949,-1.3599734,-1.32664 
439, . 4489942 5, 2. 4467B 25, -2. 91 68448E-4, 3. 1868334E- 3,-6.69886646-3, -2 
5. 453835E-2, 3. 44934686-2, .5086266,-1.9280242,2.9076944,1.6993272,-3 
6. 646D55E-5, 4. 5526191 £-4, -1.1164777E-3, -4. 907676-3,8.6233676-3, .169 
75422, -. 9640 121, 2. 907 69 44, 1.6 993272, -3.2409387E -5, 3. 983 54 186-4,-9. 5 
8698091 E-4,-4.089725E-3, 6. 8986936E-3,. 127 15665, -.64267473, 1.4538472 
9,1.6993272/ 
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DATA H I *S 1/1.176126027040027, 135*0144920564057,27.74745767553619,2 

15*21635718674784/ 

K = 1 

1 J = 1 

S=T/ 1. 003 

IF < S.GT.5.3405286D-1 > J=2 
CP=A(1,K,J> 

DO 2 N=2 * 8 

2 CP = CP*S+A(N»K»J> 

GO TO (3,4,5) ,K 

3 CP=CP*S+A<9, 1, J) 

RETURN 

ENTRY CS(T> 

K -2 

GO TO 1 

4 CP=CP*S+A(9,2,J)*DL0G(S)+S1 ( J) 

RETURN 

ENTRY CH(T) 

K=3 

GO TO 1 

5 CP=T*( CP*S*A(9,3,J ) ) «-HI ( J) 

RETURN 

END 


C DECK HOLOG 

LOGICAL FUNCTION LGFN <P,T,J,Z> 

DIMENSION Z(6) 

J-l 

LGFN=. TRUE. 

IF <P.LT.0.1.0R.P.GT.l.lE7.0R.T.LT.179.0.0R.T.GT.3010.0.0R.Z(l>.LE 

1. 0.0. DR. Z (2 ) .LE.O.O.OR.Z(3).LE.O.O) RETURN 
J=0 

LGFN=. FALSE. 

RETURN 

END 


\ 


C DECK HODATA 

BLOCK DATA 
COMMON /LDATA/ R<6> 

DATA R/l. 0,259. 8347, 31.9988, 1 . 804444E-3, 5. 6, 1 . 4/ 
END 
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C DECK HOTLG 


SUBROUTINE TLOGIC <P,T> 
IF <T.LT. 179.0) T=179.0 
RETURN 
ENO 


Normal Hydrogen 


For normal hydrogen, the equations for the compressibility factor and the ideal-gas 
specific heat are the result of curve fits of the data in reference 5. These equations 
cover a temperature range from 70 to 600 K. The pressure range extends to 100x10^ 
pascals. 

In order to check the accuracy of these curve fits, the values of Z and C calcu- 
lated by these routines were compared with those tabulated in reference 5. Over most 
of the temperature and pressure range, Z agrees to within 0. 1 percent and C p agrees 
to within 1 percent. The maximum disagreement is 0. 2 percent for Z and 5 percent 
for C p . 

The card listing for these routines follows. The prefix for the deck names is NH. 

C DECK NHZETA 


SUBROUTINE ZETA (KK,PP,TT,Z) 

DIMENSION Z ( 6 ) » BCD(6,4,2), A(4,2) 

DATA BCD/2.8086554E-3,-1.5620062E-2,3.377657E-2,-3.6316204E-2,5.87 
14617E-3,8.1391454E-3»-l. 1234621E-2? 4. 68601 85E-2» -6. 75531 39E-2 , 3. 63 
21620 4E -2,0., 8. 1391454E-3 ,-l .4043277E-2, 6 . 2480247E-2, -. 10132971 , 7. 2 
3632 409E-2 , -5 . 87461 7E-3, 0.,5.6173107E-2,-. 18744074, . 20265942 ,-7.263 
42409E-2,2*0.,-3.4569934E-5,1.5821139E-4,-2.2152582E-4,5.7982627E-5 
5,1.5996837E-4,4.0524086E-5,1.3827974E-4,-4.7463417E-4,4.4305163E-4 
6,-5.79B2627E-5,0.,4.0524086E-5,8.6424835E-5,-3.1642278E-4,3.322887 
72E-4,-5.79B2627E-5,-7,9984166E-5,0.,-3.4569934£-4,9.4926833E-4,-6. 
86457745E-4»5.7982627E-5» 2*0./ 

S=100. 3/TT 
P = PP 
K=KK 

DO 3 J-1,4 

IF ( K. EQ. 2. AND. J. NE. 2 ) GO TO 3 
DO 2 1=1,2 
A| J,I)*BCD(1,J,I) 

DO 1 M = 2 * 6 

1 AU,I|*AI J,I)*S + BCOf M,J, I) 

2 CONTINUE 

IF ( K. EQ. 1. AND. J.EO. 1 ) GO TO 4 

3 CONTINUE 

Z<2)=1.0+< AC 2, 1)+A<2,2)*P)*P 
Z(4)=( A<2,1 )+A< 2, 2)*P/2.0)*P 
IF ( K. EQ. 21 RETURN 
Z(5)=(A(3,1 )+A(3,2)+P)*P 
Z(6)*( A(4,1)+A(4,2)*P)*P 
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4 Zm«1.0HAU,i)+A(l,2)*PI*P 

Z (3) =l.0+I2.0*A{ 1,1) +3.0*A< 1,2) *P)*P 

RETURN 

ENO 


C DECK NHTEMP 

DOUBLE PRECISION FUNCTION CP (T) 

DIMENSION A ( 8, 3, 2 ) 

DOUBLE PRECISION CS, CH, S ,H I ( 2 ) , SI< 2 > 

DATA A/4.4l52734E-3,-5.1779658£-2, .1990402 ,-9.97 81 436E-2,- 1.2 598 36 
13*3. 271702,-2. 40 18526, 2. 05 33733*6. 3075334E-4»-8.629943E-3»3. 98 08 04 

2E-2*-2.4945359E-2»-. 41994543, l .635851 *-2. 401 8526 * 2.0533733* 5 . 51909 
31BE-4,-7.397094E-3,3.3173367E-2,-1.9956287E-2,-. 31495908,1.0905673 
4,-1.2039263, 2.0533733, 1.47 19026E-6, -2. 5457897E-5,-!. 6392986E-4, 4. 5 
5782 095 E-3, -1.8045494 E-2, - 6 . 6869337 E-2, .5687798,1 .5362018,2. 1027L8E 
6-7, -4. 2429828E-6*-3. 278 5972 E -5, 1.1445524E-3,-6.015l647E-3,-3.34346 

769E-2*. 5687798, 1.5362018, 1.8398783E-7, -3.636 84 24E-6, -2. 7321643E-5, 
89. 1564192E-4, -4. 51 13735E-3, -2. 2289779E-2,. 2843899, 1.5362018/ 

DATA HI, SI/ 118. 5090360344103, 85. 18245132337324, 12. 03503970042298, 1 
10.19779192852219/ 

K = 1 

1 J = 1 

S = T/ 1 » 0D2 

IF ( S.GT.3.30230786D0) J=2 
CP=All,K, J> 

DO 2 N=2 , 7 

2 CP=CP*S+A(N,K,J) 

GO TO < 3, 4, 5 ) ,K 

3 CP=CP*S+A(8,1,J) 

RETURN 

ENTRY CS(T> 

K = 2 

GO TO l 

4 CP=CP*S + A(8»2,J> *DLOG( S ) +S I ( J ) 

RETURN 

ENTRY CHIT) 

K = 3 

GO TO l 

5 CP=T*(CP*S+A<8,3,J))+HI( J) 

RETURN 

END 


C OECK NHLOG 

LOGICAL FUNCTION LGFN <P,T,J,Z) 

DIMENSION Z ( 6 ) 

J=1 

LGFN=. TRUE. 

IF ( P.LT.O. 1.0R.P.GT.1.1£7.DR.T.LT.69.0.0R.T.GT.610.0.0R.Z( D.LE.O 
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l.O.OR.Z(2).LE.O.O.OR.Z(3).LE.O.O) RETURN 
J=0 

LGFN=. FALSE. 

RETURN 

END 


C OECK NHDATA 

BLOCK DATA 
COMMON /LDATA/ R { 6 > 

DATA R/i. 0,4124. 33, 2.G1594, 1.136808E-4, 5. 6, 1.4/ 
END 


C DECK NHTLG 

SUBROUTINE TLOGIC IP,T) 
IF (T.LT.69.0) T =69. 0 
RETURN 
END 


Parahydrogen 


For parahydrogen, the equation for the compressibility factor is from reference 9. 
For temperatures below 100 K, this equation is based on the parahydrogen data of refer- 
ence 10. For temperatures above 100 K, this equation is based on normal-hydrogen data 
such as that found in reference 5. The equation for the ideal-gas specific heat is the 
result of a curve fit of the data in reference 11. For temperatures from 13 to 100 K, the 
pressure range extends to 300x10^ pascals. For temperatures from 100 to 400 K, the 
pressure range extends to lOOxlO 5 pascals. 


Except for the critical region, reference 9 estimates the mean error in Z to be 
0. 3 percent. Over most of the range from 100 to 400 K, the values of Z calculated by 
these routines agree with those tabulated in reference 5 to within 0. 1 percent. No com- 
parison of C was made. 

The card listing of these routines follows. The prefix for the deck names is PH. 
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C DECK PHZETA 


SUBROUTINE ZETA( KK,PP,TT,Z ) 

DIMENSION Z ( 6 ) 

DATA A,B1 ,B2,B3,B4,B5,B6,B7, B8 , B9, BIO, B1 1 , B12 , B1 3 , B14, Bl 5/-4. 4291 0 
1289E-4,1.0229558E-2,-7.81250336E-l,-43.7354838,700.54484,-61088.73 
208, 7. 8151 9168E-5, 7. 73767544E-3, 3. 58441 339E-7, -32. 0855334, 3046. 6492 
34,-581 06.0523, 2 .84494 81 6E- 3, -4.98394981E-1, 10.78378 94, 1.9225229E-8 
4/ 

C0(X,Y,Z)=X*S3+Y*S4+Z*S5 

K=KK 

P=PP 

T=TT 

Sl=1.0/T 

S2=S1*S1 

S3=S1*S2 

S4=S1*S3 

S5=S1*S4 

P2=P*P 

EXPO=EXP< A*PZ) 

AA=1.3Z(2.0*A> 

BB=l.3/A 

Cl=CO( B9, BIO , Bl 1 ) 

C2=CO( 812»B13»B14) 

IF (K.EQ.2) GO TO 1 
BA1=B1+B2*S1+B3*S2+B4*S3+B5*S5 
B A2 = B6+B7*S 1 

ZA=1 .0+( BA1+(BA2+IB8+B15*S1*P2)*P)*P)*P 
ZB=EXP3*P2*(C1+C2*P2 ) 

ZU)= ZA+ZB 

ZA=1.0+(2.0*BAl+(3.0*BA2+( 4. 0+B8+6. 0*B15*S1*P2 ) *P ) *P > *P 
ZB=P2*EXPO*(3.3*Ci+< 2.0*A*C1+5.0*C2+2.0*C2*P2*A)*P2) 

Z ( 3 ) = ZA+ZB 
IF (K.EO.l) RETURN 

1 BT=B1-B3*S2-2.0*B4*S3-4.0+B5*S5 

C1P=-CD(3.0*B9, 4.0*B10,5.0*B11) 

C2P=-CD( 3.0*812, 4. 0*B1 3, 5.0*B14) 

CC1=C1+C1P 

CC2=C2+C2P 

ZA s i.0+(BT+(B6+88*P)*P)*P 

Z8=P2*EXP0*(CC1+CC2*P2) 

Z C 2 > — ZA+ZB 

ZA=(BT+(B6/2.Q+B8*P/3.0I*P)*P 
ZB=AA+(EXPO+(CCl + ( P2-BB)*CC2 >+B8*CC2-CCl) 

Z(4> = ZA+ZB 

IF (K.EQ.2) RETURN 

C1PP=CO(9.D*B9,16.0*B10,25.Q*B11) 

C2PP=CO(9.0*B12, 16. 0*B13, 25.0*B14) 

ZA=-((B2*$1+2.0*B3*S2+3.0*B4*S3*5.0*B5*S5)+(B7*S1/2.0+B15*S1*P**3/ 
15.0) *P)*P 

ZB=AA*(EXP0*(C1P+IP2-BB) *C2P)+BB*C2P-C1P) 

Z 1 5 ) = ZA+ZB 

CC1P=C1P+C1PP 

CC2P=C2P+C2PP 

ZA=(2.3*B3*S2+6.0*B4*S3+20.Q*B5*S5)*P 

ZB=AA*(EXP0*(CC1P+(P2-BB)*CC2P)+CC2P*BB-CC1P) 

Z ( 6 I = ZA+ZB 

RETURN 

END 
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C OECK PHTEMP 


DOUBLE PRECISION FUNCTION CP<T> 

DIMENSION A<6,3,4) 

DOUBLE PRECISION SL< 3 ) , S 1 1 4 » ,HII 4 ) , S 

DATA A/5*0. 0, 1.5, 5*0. 0,1. 5, 5*0.0, 1.5, D.Q ,-5.477243, 14. 151817, -10.9 

192694. 3. 4450167. 1.118 1079. 0. 0, -1.36931 075, 4. 71 727233, -5. 496347, 3. 4 

2450167.1.1181079.0. 0,-1.0954486,3.53795425,-3.66423133,1.72250835, 
31. 11 81 079, -.9606211, 8. 1602445, -26. 350196, 39. 1680848, -25. 0299192, 7. 
42610 8881, -.19212422, 2.04006112,-8.78339867, 19.5840424,-2 5.0299192, 
57. 26108881, -.1601035 17, 1.6320489,-6.587549, 13.0560283, -12.5149596 
6, 7. 261 08881, 2*0. 0,-2. 3415379E-2, .3216503,-1.4584629,4.7141156,2*0. 
70,-7.80512633E-3,.16082515,-l.4584629,4.7141156,2*0.0,-5.85384475E 
8-3,. 107216767, -.7292 3145, 4. 7141 156/ 

DATA S I, H 1/9.9827740 5790 1375, 8. 867877636855692,2 2. 54596376458948, 
11 1.1 4950174192363, -.02542447420268701, 3. 222638996244185, -103. 532 18 
278709574,-282.1978829444169/ 

DATA SL/. 299951 327D0, .8020822700, 2. 0023983100/ 

K = 1 

1 S=T/1. 0D2 

I F ( S.GT.SLI 1 ) } GO TO 2 
N = 1 

GO TO 5 

2 IF(S.GT.SL(2MGO TO 3 
N=2 

GO TO 5 

3 IFIS.GT.SLOMGO TO 4 
N=3 

GO TO 5 

4 N=4 

5 C P= A ( 1 ,K, N) 

DO 6 1=2,5 

6 CP=CP*S+A(I ,K,NI 
GO TO ( 7, 8, 9 ) , K 

7 CP=CP*S+A(6»1,N) 

RETURN 

ENTRY C S C T 1 
K = 2 

GO TO 1 

8 CP=CP*S*A(6,2,N)*DLOG(S)+SI < N ) 

RETURN 

ENTRY CHIT) 

K = 3 

GO TO 1 

9 CP=T*(CP*S+A(6»3»N>l*HI(N) 

RETURN 

END 


C DECK PHLOG 


LOGICAL FUNCTION LGFN (P,T,J,Z) 
COMMON /LDATA/ XKV,R, XMW ,RC, D2, G 
DIMENSION Z ( 6) , A<8) 
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S*T/IOD.O 
J = 1 

LGFN=. TRUE. 

IF (P.LT.0.1.0R.P.GT.3.1E7.0R.S.GT.4.1.0R.S.LT.0.13.0R.Zm.LE.Q.O 
1.0R,Z(2>.LE.O.O.OR.Z(3I.L£.O.O.OR.<P.GT.l.lE7.AND.S.GT.l.ll J RETUR 

2N 

IF ( S.GT. 0.329761 GO TO 2 
PL=A(1) 

DO 1 N=2 1 7 

1 PL=Pl*$+A(N ) 

PL=XKV*EXP( PL+AI8 ) / S ) 

IF (P.GT.PL) RETURN 

2 J-0 

LGFN=. FALSE. 

RETURN 

END 


C DECK PHDATA 

BLOCK DATA 
COMMON /LDATA/ R( 6) 

DATA Ra.O,4124.33,2.01594,1.136808E-4, 5. 6, 1.4/ 
END 


C DECK PHTLG 

SUBROUTINE TLOGIC<P,T> 

DIMENSION A ( 6 ) 

DATA A/5.6592537E-5,-3.5285986E-4,-5.8183381E-3,9.6000865E-2»6.523 
1328BE-1, 3. 618144/ 

IF(T.GT.32.976)RETURN 
V=ALOG( PI 
S=A ( 1 ) 

DO 1 N=2?6 
I S=S*V*A(N) 

IF(T.LT*S)T=S 

IF(T.LT.13.0)T=13.0 

RETURN 

END 
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Helium 


For helium, the equation for the compressibility factor is from reference 12. Since 
helium is a monatomic gas, the value of the ideal-gas specific heat at constant volume 
C y jcj/R i s 1-5. Reference 12 claims a temperature range from 3 to 300 K and a pres- 
sure range extending to 101x10^ pascals. This report assumes a temperature range from 
5.4 to 400 K and a pressure range extending to 300x10^ pascals. These are the same 
ranges assumed in reference 6, which uses the same equations. 

For pressures to 101x10^ pascals and temperatures from 10 to 300 K, reference 12 
estimates that the calculated values of p, H, and S are accurate to within 3 percent and 
that the values of C p are accurate to within 5 percent. 

The card listing of these routines follows. The prefix for the deck names is HE. 

C DECK HEZETA 


SUBROUTINE ZET A { KK, PP, TT , Z ) 

DIMENSION Z<6) 

DATA A'Bl,B2,B3tB4,B5»B6'B7 v B8.B9.B10.Bll,B12»B13,B14,Bi5 > BL6/ 
l-4.057E-4,4.0665013E-3,-l. 1267764E-1. 2. 3039266E- 2 , -5. 74688186-2 , 
21.3691 368E-1, 9. 7390626E-6, 7.0543876E-4, -5. 3854984E-6.-3. 8053762E-3 
3,2.625l79E-2,7.674266E-2,-8.790491lE-7, 1. 9960611E-6, -8. 1 300167E-6, 
43.67435B3E-8,-3.4049435E-ll/ 

CO(X,Y r Z)=X*S3+Y*S4+Z*S5 
K=KK 
P = PP 
T=TT 

S1=1.0ZT 

S2=S1*S1 

S3=S2*S1 

S4=S3*S1 

S5=$4*S1 

AA=A*Sl 

P2=P*P 

EXPQ-EXP(AA*P2) 

C1=CO(89»B10*B11 > 

C2=C0( B12,B13,B14) 

IF (K.E0.2) GO TO l 

BA1 = B1«-B2*S1 + B3*S2+B4*S3+B5*S5 

BA2=B6+B7*S1 

ZA=1.0+(RAI*(BA2+(BB*S1+(B15*S1+B16*S1*P)*P)*P)*P )*P 

ZB=ICH-C2*P2)*P2*EXP0 

ZC1»* ZA+ZB 

ZA=1.0M2.0*BAl + ( 3.0*BA2 + (4.0*B8*Sl + < 5.0*B15*S1 + 6.0*816*S1*P)*P)*P 
1>*P>*P 

ZB= I 3. 0*C1+ I 2.0*AA*C 1+5.0*C2+2.0*AA*C2*P2)*P2)*P2*EXP0 
Z ( 3 ) = ZA+ZB 
IF (K.EQ.l) RETURN 

1 BT=Bl-B3*S2-2.Q*B4*S3-4.0*B5*S5 

C1P=-COI3.0*B9,4.0*B10, 5.0*811 > 

C2P=-C0I 3.0*812 » 4. 0*B13, 5.0*B14> 

Cll=(2.0*Cl+CiP)/AA 

C22=(3.0*C2+C2P)/AA 

ZA=1.0*(BT+B6*P)*P 

ZB = P2*EXP0*IC1*C1P+I -AA*C1+C2+C2P-AA*C2*P2 »*P2> 
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Z 1 2 ) = ZA+ZB 
ZA=IBT+B6*P/2.0)*P 

ZB=D.5*(EXP0*(Cll-C22/AA+( -C 1+C22-C2*P2 > *P 2 >-Cll+C22/AA ) 

Z I A ) * ZA+ZB 
IF (K.EQ.2) RETURN 
Cll = ( CI+C1P ) f AA 
C22=I2.0*C2+C2P)/AA 

Z A= ( -( B2*S1 +2.0*83*S2+3.0*B4*$3+5. 0*B5*$5 ) ♦( -B7+S1/2. 0+(-B8*Sl/3.0 
l+(-B15*Sl/A.D-B16*Sl*P/5.0)*P)*P)*P)*P 
ZB=0.5*IEXPO*ICll-C22/AA+< -C 1+C22-C2*P2 )*P2 )-Cil+C22/AA) 

H5)~ ZA+ZB 

CIPP=COI9.0*B9,16.0*B10,25.0*BU> 

C2PP=C0I 9.0*812, 16. 0*B13 ,25.0*814) 

C11=(2.0*C1+3.0*C1P+C1PP )/ AA 
C22=I6.0*C2+5.0*C2P+C2PP )/AA 
ZA= ( 2.0*B3*S2+6.0*B4*$3+20.*B5*S5)*P 

ZB-0.5*(EXPO*(Cll-C22/AA+l -2.0*1 C1+C1P ) +C22+ I AA*C1-3.0*C2-2.D*C2P+ 
IAA*C2*P2)*P2)*P2 )-Cl 1+C22/AA) 

Z(6) = ZA+ZB 

RETURN 

END 


C DECK HETEMP 

DOUBLE PRECISION FUNCTION CP<T) 

DOUBLE PRECISION $,C15,SI,HI 

DATA C15»SI »HI/1.5D0, 4. 7506300,6. 98973DQ / 
CP = Cl 5 
RETURN 
ENTRY CSI T ) 

S = T 

CP= CI5*OLOG( S ) +S I 

RETURN 

ENTRY CHIT) 

S = T 

CP= Cl 5*S+HI 

RETURN 

END 


C DECK HELOG 

LOGICAL FUNCTION LGFNIP, T,J,Z) 

DIMENSION Z ( 6 ) 

J = l 

LGFN= .TRUE. 

IF (P.LT.0.1.0R.P.GT.3.1E7.0R.T.LT.5.4.0R.T.GT.410.0.0R.Zm.LE.O. 
IO.OR.ZI2I.LE.O.O.OR.ZI3) .LE.O.O) RETURN 
J=0 

LGFN= .FALSE. 

RETURN 

END 



C DECK HEDATA 


BLOCK DATA 
COMMON /LDATA/ RI6) 

DATA R/l. 0,2077. 25* 4. 0026,2.5386636-4,. 3333,1. 6666667/ 
ENO 


C DECK HETLG 

SUBROUTINE TLOGIC(P,T) 
IFIT.LT.5.4) T = 5.4 
RETURN 
END 


Argon 


For argon, the equation for the compressibility factor is from reference 13. Since 
argon is monatomic, the value of the ideal-gas specific heat at constant volume C v.i</ R 
is 1.5. Reference 13 claims a temperature range from 86 to 300 K and a pressure range 
extending to 1010x10^ pascals. This report assumes a temperature range from 80 to 
400 K and a pressure range extending to 1000x10^ pascals. 

The values of Z and calculated by these routines (using the state equation of 
ref. 13) were compared with those tabulated in reference 5. Over most of the tempera- 
ture range 200 to 400 K and for pressures under 101X10 5 pascals, Z agrees to within 0. 1 
percent and agrees to within 1 percent. The maximum disagreement (at about 200 K 
and 100 bar) was 2. 4 percent in Z and 15 percent in Cp. 

The card listing of these routines follows. The prefix for the deck names is AR. 

C DECK ARZETA 

SUBROUTINE ZETA( KK,PP,TT ,Z ) 

DIMENSION Z<6» 

DATA A,Bl,82,B3,B4,B5,B6»B7,B8,B9,B10,Bll,B12,Bi3,&14,B15/-3.0D798 
1 3 66 E -6 ,7. 92 5 59677E-4, -2* 73770 136E- 1,-20. 5241378,-808.296022 .297857 
26.6*5. 382625E-7,-3.57177204E-5, 4.328578956-10,3.67665496*493. 10182 
35 *-87715.9942,-3.103016966-6, 2. 22667069E-3, 5. 232 79441E-2, 8. 3194622 

47E-14Z 

C0( X,Y,Z)=X*S3+Y*S4+Z*S5 
K=KK 
P = PP 
T*TT 

Sl=l.0/T 

S2=S1*S1 

S3=S1*S2 

S4=S1*S3 
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S5=S1*S4 

P2-P*P 

EXPO=EXP< A*P2> 

AA=1.0/(2.0*A) 

BB=l »0/A 

C1=C0( B9,B10,Bll) 

C2=C0<B12,B13»B14) 

IF (K.EQ.2) GO TO 1 

BAl=8i+B2*Sl+B3*S2+B4*$3+B5*S5 

BA2=B6+B7*S1 

ZA=1.0+<BAl+(BA2+(B8+Bl5*Sl*P2)*P)*P»*P 

ZB=EXP0*P2*(Cl+C2*P2> 

Z(D= ZA+ZB 

ZA=1.0+(2.0*BAl+< 3.0*BA2+( 4. 0*B8+6.0*B15*S1*P2 )*P ) *P > *P 
ZB=P2*EXPO*< 3.0*C1 + (2.0*A*C1 + 5.0*C2+2.0*C2*P2*A)*P2> 

Z ( 3 J = ZA+ZB 
IF (K.EQ.D RETURN 

I BT=B1-B3*S2-2.0*84*S3-4.0*85*S5 

C1P--C0( 3. 0*89, 4.0+B10, 5.0*B11 » 

C2P=-CO(3.3*812,4.0*B13,5.0*B14) 

CCI^CI+CIP 

rr?=r2+r2P 

ZA=1.0+(BT+(B6+B8*P)*P»*P 

Z8=P2*EXP0*(CC1+CC2*P2» 

Z(2>* ZA+ZB 

ZA=(BT+<B6/2.0+B8*P/3.0)*P)*P 
ZB=AA*(EXPO*(CCl+(P2-BB>*CC2) +BB*CC2-CC1 ) 

Z<4) = ZA+ZB 

IF ( K. EQ. 2 ) RETURN 

C I PP=C 019,0*89, 16. 0*610, 25.0*Bll) 

C2PP=COI9.0*B12,16.0*B13,25.0*B14» 

ZA--nB2*Sl+2.0*B3*S2+3.0*B4*S3+5.0*B5*S5) + IB7*Sl/2.0+B15*Sl*P**3/ 

15,D)*PI*P 

ZB=AA*(EXP0*<CIP+(P2-BB»*C2P)+BB*C2P-C1P) 

Z ( 5 > - ZA+ZB 

CCIP=CIP+CIPP 

CC2P=C2P+C2PP 

Z A= ( 2. 0*B3*S2+6.0*B4*S3+20. 0*B5*S5 ) *P 
ZB=AA*1EXP0*(CC1P+(P2~B8)*CC2P )+CC2P*BB-CClP » 

Z<6)= ZA+ZB 

RETURN 

END 


C DECK ARTEMP 

DOUBLE PRECISION FUNCTION CPIT) 

DOUBLE PRECISION $,C15»SI 

DATA C15.SI/1. 500, 10. 54223318/ 

CP=C15 

RETURN 

ENTRY CS(T) 

CP=C15*DL0G( DBLE(T) ) +SI 
RETURN 
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ENTRY C-H ( T ) 
CP-C15*T 
RETURN 
END 


C DECK ARLOG 

LOGICAL FUNCTION LGFN (P,T,J t Z) 

COMMON /LDAT A/ XKV , R , XMW ,RC» D2 , G 
DIMENSION Z<6) 

S=T/103.0 

J=1 

LGFN=* TRUE. 

IF (P.LT.0.1.0R.P.GT.1.1E8.0R.S.GT.4.1.0R.S.LT.0.8.0R.Z< D.LE.0.0. 
10R.Z(2).LE«O.O.OR.Z(3).LE.O.O) RETURN 
IF 1S.GT. 1.5086) GO TO 1 

PL=XKV*EXPl-9. 3287334/ S+ 24. 471369+1 -3. 5108148+1. 0598825*$ ) *S ) 

IF (P.ST.PL) RETURN 
i J = 0 

LGFN=. FALSE . 

RETURN 

END 


C DECK ARDATA 

BLOCK DATA 
COMMON /LDATA/ R<6> 

OAT A R/l.0,208.1306, 39. 948,2. 533716E-3, 5. 6, 1.6666667/ 
END 


C DECK ARTLG 

SUBROUTINE TLOGIC(P»T) 

DIMENSION ACS) 

12535 E_3 »“l* 3007 3 85E-1 , 3.9880245E-1 ,12.49659,-24.58815/ 
IFIT.GT. 150,86) RE TURN 
V=ALGG( P) 

S = A ( 1 ) 

DO 1 N=2, 5 
1 S=S*V+A( N ) 

IF(T.LT.S)T=S 
IF< T.LT.80. )T=80. 

RETURN 

END 
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Steam 


For steam, the equations for the compressibility factor and the ideal-gas specific 
heat are derived from equations (1) to (3) in reference 14. The temperature range for 

5 

these equations is from 273 to 1600 K. The pressure range extends to 1000x10 pascals. 

It is estimated, that except for the critical region, the values of Z are accurate to 
within 0. 2 percent. 

The reference state for the internal energy and entropy is that of the liquid at the 
triple point. (The value of the enthalpy at this state is 0. 6 joules per kilogram.) 

The card listing of these routines follows. The prefix for the deck names is ST. 

In this case, Z(l) to Z(6) are calculated in subroutine ZETA regardless of the value 
of K. 

C DECK STZETA 

SUBROUTINE ZETA (K,PP,TT,Z> 

DIMENSION Z(6), A(8,7), 0(7, 7), PJU0,2), PAJ(2), FI(7,3> 

DOUBLE PRECISION D,PD 

DATA A/29.492937,-132.13917, 274.64632,-360.93828,342.18431,-244.50 
1042,155.18535,5.9728487,-5.1985860,7.7779182,-33.301902,-16.254622 
2,-177.31074,127.48742,137.46153,155.97836,6.8335354,-26.149751,65. 
3326396, -26.1 81978, 4*0.0,-. 15641040, -.72 546108, -9. 2734289, 4. 3125840 
4,4*0.0,-6.3972405,26.409282, -47 .740374, 56. 3231 30, 4*0. 0,-3. 96 61401, 

51 5. 453061, -29. 142470, 29. 568796, 4*0.0, -.69048554, 2. 740741 6, -5. 102 80 
670, 3. 9636085, 4*D.0/ 

DATA D/-410. 3084800, 337. 31 18000, -137. 46618D0, 6. 7874983D0, 136. 87317 
1D0, 79. 84797000, 13. 04 1253D0, -416. 0586000, -209. 88866D0, -733. 968480 0, 
210.40171700,645.8188000,399.175700,71.53135300,1137.363504,-2038.8 
373960,-808.0992960,-11.77655784,634.6463840,415.0811440,80.4646916 
40,1997.081280, 1007.465568,3523.048704,-49.92824160,-3099.930240,-1 
5916.043360,-343.3504944, 3106 .844208,-3657. 9706D0 ,-148. 2616320 ,-44. 
635654958, -22. 344832D0, 31. 81088800, 17.86667720,-1465.1822592,11801. 
7526144,10924.9740288,-43.329005568,-9246.1631232,-5824.4762112,-10 
872. 931 50848, -9585. 990 144, -4835. 8347264, -16910. 6337792, 2 39. 65555968 
9,14879.665152,9197.008128, 1648.08237312/ 

DATA TC,E«PAJ, PJ/ 1. 54491 20 ,-4.8, .634, 1.0,20*1.0/ 

P=PP/1.0E3 
T=1 . 0E3/TT 
DO 1 N=1 , 2 
PJ(4,N»=P-PAJ(N) 

PJ(2,N)=1.0/PJ(4,N) 

PJ(l,N)=PJ(2,N)**2 
DO 1 1=5,10 

1 PJ( I ,N)=PJ( I— 1»N)*PJ ( 4»N) 

EXPO=EXP(P*E> 

PD=DBL E ( P ) 

N = 1 

00 3 J=1 ,7 

F I ( J,i)=EXPO*(SNGUD( J,1 > + D< J,2)*POI 1 

FI ( J,2)=EXP0*( SNGUDI J, 1 >♦! DC J, 3)*D( J,4)*PD»*PD) I 

FI < J,3)=EXPO*(SNGL(D( J,5)+(0< J,6>+0< J,7)*PD)*PD > ) 

DO 2 1=1,8 

IF ( J.GT.2.AND. I.GT.4) GO TO 3 
FI ( J,i)=FI( J, 1 >+A( I, J)*PJ< 1 + 2, N ) 

X I = I 


67 



FI< J»2)=FI( J,2)+AU,J)*(XI*P-PAJ(N> >*PJI I+i,N) 

2 FI(J,3) = FI(J,3H-A<I,J)*( Xl-1.0)*UI*P-2.0*PAJ< N> )*PJ(I»N) 

3 N=2 
PT=P*T 
SC=T-TC 

Q=FI(l,lJ+Fl(2,l»*SC 

PQ=FN1,2)+FI(2,2)*SC 

PPQ=FI(l»3)+Fl(2»3)*SC 

VV=SC 

V = 1.0 

S=T-2.5 

$2=S*S 

T Q~F 1(2*1) 

TTQ=0. 0 
PTQ=FI (2*2) 

DO 4 J-3,7 

V=V*S 

VV=VV*S 

Q=Q+W*FI ( J,l) 

PQ=PQ*VV*FI< J,2) 

PPQ= PPQ*VV*F I ( J, 3 ) 

X J= J-2 

Tl=( XJ*SC+S >*V/S 
TQ=TQ+T1*FI ( J,l) 

PTQ = PTQ«-T1*FI(J»2) 

4 TTQ = TTQ-t-XJ*( ( X J-1.0 ) *SC+2.0*S )*V*Fll J, 1 )/S2 
Zm=l.O+P*PQ 

Z(3)=2.0*Z(1 )-1.0+P*P*PPQ 

Z(2)=Z(1)-PT*PTQ 

Z(5)=-PT*TQ 

Z ( 4 ) =P*Q+ Z ( 5 J 

Z (6 ) =T +PT *TT Q 

RETURN 

END 


C OECK STTEHP 

DOUBLE PRECISION FUNCTION CP(T) 

DIMENSION A ( 6 f 3 ) 

DOUBLE PRECISION S 

DATA A/-. 21 D28060,. 5 3437455, -.47667309, 1.8177938, 2. 1911746,. 099672 
1813,-. 052 5701 50,. 178 12485, -.23833655, 1.8177938, 2. 191 1746, -.0996728 
21 3, -.0 42056120, . 13359364, -.15889103, .90889688, 2.1911746,99.672813/ 
DATA HI, SI/4612. 7316, 17.238170/ 

K = 1 

1 S=T/ 1. 0E3 
CP- A ( 1 ,K> 

DO 2 1=2,4 

2 CP=CP*S+A ( I ,K ) 

GO TO (3,4,5) ,K 

3 CP=CP*S+A<5,K)+AI6,K)/S 
RETURN 

ENTRY CS( T ) 

K=2 

GO TO 1 
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4 


CP=CP*S+A(5,2)*DLOG< S ) +A ( 6, 2 I /S+S I 
RETURN 
ENTRY CHIT) 

K=3 

GO TO 1 

5 CP = T*(CP*S+A(5,3> )+A ( 6, 3 )*DLOGI S ) +HI 

RETURN 
ENO 


C DECK STLOG 

LOGICAL FUNCTION LGFN (P,T,J,Z) 

DIMENSION Z ( 6 ) • A<9) 

COMMON /LOATA/ XKV, R , XMW ,RC, D2, G 

DATA A/2.828160iE-5,-5.8328303E-4,4.084426E-3,-5.9629737E-3,-7.000 
13883E-2,5.6897698E-1 ,-3.091692,35.493889,-64.278808/ 

S=T/10D.O 

J=1 

LGFN=. TRUE. 

I F <P.GT.1.01E8.0R.P.LT.0.1.0R.S.GT.16.1.0R.S.LT.2.73.0R.Z( D.LE.O 
1.0.0R.ZI2 I.LE.O.O.OR.Zi 3KLE.O.0) RETURN 
IF (S.GT. 6. 472861 GO TO 2 
PL=A ( 1 ) 

DO 1 N=2 , 8 

1 PL=PL*S+A(N I 

IF (P.GT.XKV*EXP(Pl«-AI9>/S) > RETURN 

2 J = 0 

LGFN=. FALSE. 

RETURN 

END 


C DECK STOAT A 

BLOCK DATA 
COMMON /LDATA/ R(6) 

DATA ft/1. 0,461.51, 18. 01534, 9. 82042E-4, 5. 6, 1.3333333/ 
END 


C DECK STTLG 

SUBROUTINE TLOGIC (P,T) 

DATA Al, A2, A3,A4, A5, A6/2 .6832873, -0. 1739764, . 038 877378 ,-l . 72571 2 7E 
1-3, -9. 1550282E-6, 3.4056076E-6/ 

S=T/1D3.0 

IF I S. GT. 6. 47286 ) RETURN 
V=ALOG(P) 

SS = A1M A2+{ A3+( A4+< A5+A6*V I *V ) *V I *V ) *V 



Methane 


IF (S.LE.SS) S=SS 

IF (S.LT.2.73) S=2.73 

T=100.5*S 

RETURN 

ENO 


The routines for methane are described and presented in reference 3. The compres- 
compressibility -factor equation is from reference 15, and the ideal-gas specific -heat 
equation is the result of a curve fit of the data in reference 16. The temperature range is 
from 70 to 600 K and the pressure range extends to 400x10 pascals. 

Even though the card listing of the methane routines is in reference 3, this listing 
is included in this report for convenience. The prefix for the deck names is ME. The 
subroutine ZETA is not independent, but the combination of subroutines ZETA and 
POLY is independent. 

C DECK MEZETA 

SUBROUTINE ZETA <KK,PP,TT,Z) 

COMMON /VALUE/ F( 4, 4 ) » G( 6, 4) 

DIMENSION Z ( 6 ) 

DOUBLE PRECISION F , G , B1 , B2 , B3, B4, B5 , A1 , A2, A3, A4, A 5 , El , E2 , P A , T A , THl 
l,TH2,TH3,TH4,0l,D2,D3,D4,D5,Fl,F2, UA, P, T, P l , P2 , U, T1 , RC , E XPC , R8 , E XP 
2B,ZB1, ZC1 , AB1,AB2,AB3,AB4, AB5, ZA, ZB, ZC , RBI , EXP81 , S , SS , PS 1 1 , PS 1 2, PS 
3I3,PSI4,RC1,EXPCI,PSI5,PSI6,PSI7,PSI8 
OAT A B1,B2, B3,84,B5, A1,A2, A3 , A4, A5 , E 1 , E2 , P A/4. 91 473574991686D-I 3 , 7 
1.37664223478550D-06,-l.l4587843032923D-07, 5. 8951 020951 l 1410-10, -5. 

27 43 8 22 8 13435 32D-1 3,-2. 2 3 98 31 99201 862 D00, 1. 3433 12 5 3741 270D-03, 2.759 
31 01 829 56551 0-05, -1.655469770535420-07,2. 34124562 687064D-10, -4. 60C2 
400000050000-02,-2.117700000000000-10,1.13318000000000002/ 

DATA T A, TH1, TH2,TH3,TH4, D1,02,D3, 04, D5,Fi,F2,UA/1.4771O55O0C00: 000 
12 , 1. 09934666473654D-1 4, 1 . 64873321 2840 64D07 , 1.07243639762491008,3.6 
2644888 82455 14D-1 5, -3.97760537104600000,-1. 506225 16081 086D-02, 4. 3 29 
3407407326480-04,-1.85355607372189D- 56, 2. 0528631530331 4D-09, -1.37 87 
49330005000003,1. 3441 846000 OOOODOO, 1.455 11293919343D06/ 

K = KK 
P = PP 
T = TT 
P1=P+PA 
P2=P1*P1 
U=P2*P1 
T 1 = T + T A 

RC=(Fl+F2*P)/T 
E XPC=DEXP ( RC ) 

RB= I El *-E2*U ) 

EXPB=DEXP(RB*T1> 

ZB1= (TH1*P*P2*(U-TH2 > *( TH3-U )*EXP8 ) /T 

ZC1 = ( Dl + ( D2 + (D3*<04<-D5*P ) *P ) *P ) *P ) *P*EXPC / T 

IF (K.EQ.2) GU TO 1 
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AB1=B1+A1/T 

AB2=B2+A2/T 

AB3=B3+A3/T 

AB4=84+A4/T 

AB5=B5+A5/T 

ZA = 1.3i)0+(ABl+(AB2+( AB3+ ( AB4+AB5*P ) *P )*P >*P ) *P 
Z<1 >=ZA+ZB1+ZC1 
IF (K.EQ.O) RETURN 

ZA=l.DDO+(2.0DO*A8l+ ( 3.0DO*AB2+( 4.0DO*AB3+( 5.0D3*AB4+6. ODO*A85*P )* 
1P)*P)*P)*P 

ZB=ZB1*(2.0DO*(I.ODO+P/P1) +3.0DO*P*P2*( E2*T1+1 .ODO/ < U-TH3 ) + 1 . ODOZ ( 
1U-TH2) ) > „ 

ZC=I F2*P*ZC1+EXPC*( 2 .ODO*Dl+ ( 3. 000*02 ♦( 4.000*03+ I 5.000*04+6. 00d*D5 
1*P)*P)*P)*P)*P)/T 
Z (3 )=ZA+ZB*ZC 
IF (K.EQ.l) RETURN 
I RB1=E1+E2*UA 

EXPB1=DEXP(RB1*T1 ) 

ZA=1.0D0+(Bl+<B2+(B3+( B4+B5*P)*P)*P)*P>*P 

ZB=RB*T*ZB1 

ZC=-RC*ZC1 

Z ( 2 ) =Z A+ZB+ZC 

S=E2*T 1 

SS=F2/T 

CALL POLY ( 1,1,U,T,S ) 

CALL POLY ( 1 ,2, P » T,SS ) 

ZA=MBl + IB2/2.OD0+<B3/3.0D0 + < B4/4. 0D0+B5*P/5.0D0) *P ) *P ) *P >*P 
PSIi=F(l,l)-F<2,l)+F(3,l)-F«4,l> 

PSI2=F(1,2)-F(2,2)+F(3,2)-F<4,2) 

ZB=TH4*(PSI1 *EXPB-PS I 2*EXPB1 ) 

PS13 = G(1,1)-G( 2,1)+G( 3, l )-GI 4,1)+G( 5,1)-G(6,1> 
PS14=G(1,2)-G(2,2)+G<3,2>-G(4,2)+G(5,2)-G(6,2> 

RC1 =F1 ZT 
EXPC1=DEXP(RC1) 

ZC=(PSI4*EXPC1-PSI3*EXPC)/T**2 

Z < 4) =ZA+ZB+ZC 

IF ( K. EQ. 2 ) RETURN 

CALL POLY <2,1,U,T,S> 

CALL POLY <2,2,P,T,SS) 

PSI5 = FU,3)-F<2,3)+FC3,3)-F<4,3> 

PSI6=FI 1,4) -F< 2,4)+F( 3,4)-F(4,4) 

PS17=G( 1,3>-G( 2, 3)+G( 3, 3 l-GC 4, 3 ) +G( 5, 3) -G( 6, 3 ) 

PSIB=G( 1, 4 ) — G I 2,4)+G( 3,4)-G( 4,4)+G( 5,4)-G( 6,4) 

ZA=-( Al+I A2/2.0D0+I A3/3. 0D0+ ( A4/4.000+A5*P/5.000 ) *P ) *P I *P > *P/T 
ZB-TH4*(PSI5*EXPB-PSI6*EXPBi) 

ZC=( PSI8*EXPC1-PSI7*EXPC )/T**2 
Z ( 5 I =Z A+ZB+ZC 

ZB=TH4*T*(£XPB*(RB*PSI1-(F(1,1)-2.0DO*F(2,1>+3.000*F(3,1)-4.0DO*F( 
14,1 ) )/Tl)-£XPBl*(R81*PSI2-(F(l»2 )-2.ODO*F(2,2)+3.0OO*F(3,2)-4.ODO* 
2F!4,2) )/Tl> ) 

ZC=(EXPC*< (2.0D0+RC)*PSI3-G( 1,1) +2. 3D0*G( 2 , 1 )-3. ODO*G( 3»i)+4.0D0*G 
1(4,1 1-5. *G( 5»1)+6.0D0*G(6, 1) >-EXPCl*U 2 . 000+RC1 ) *PSI 4-GC 1 , 2 ) +2.000 
2*G<2,2 >-3.0DO*G(3,2)+4.0DO*G(4,2)-5.0DO*G( 5 , 2 ) +6. ODQ*G( 6, 2 > ) >/T**2 
Z (6 ) =ZB+ZC 
RETURN 
END 
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C OECK MEPOLY 


SUBROUTINE POLY ( J»K , PP, TT, CC ) 

COMMON /VALUE/ F ( 4,4 ) ,G( 6, 4 ) 

DIMENSION A( 7 , 2 ) , B( 16,2 ) 

DOUBLE PRECISION PP, TT , CC, F, G, A, B, AA, AB, UA, D1 »D2 , D3,D4,D5,P, T,C1 , C 
12,C3,C4,C5,C6,V1,V2,V3 

DATA B/5. 4842956356422 4D03, 1 .54210957155370D01,-6.17182239272157D- 
101, 3. 13762297202746D-03, -5. 32199549075239D-06, 2. 75942703621 4580-09 
2,-1 . 23436447854431DD0, 9.41 286891 608238D- 03, 1.882 57378321 648D-02, -2 
3. 128 7931 9630D95D-05, -6. 38639458890286D- 05, -1,277 2789 177805 7D-04, 1. 
437971 351810729D-08,5.51885407242917D-08, 1.655656221728750-07,3.311 
5312443 45750D-0 7, 5*0. 0000,2. 7 59427036214580-09, 6*0. 00 00, 1.37971 351 8 
61 07290-08,5. 518854072429 17D-08, 1.655656221 72875D-07.3. 311312443457 
7500-07/ 

DATA A,AA,AB,UA,D1,D2,D3,D4, D5, VI, V2, V3/8. 13389656644895013,-5. 31 7 
142860649802006,1.979949208266470-02,2.117700000000000-10,3.9598984 
21 653293D-02 , 6. 3531000000 PO 00 0-10, 1.27 062000000000D-09, 3*0. ODOO, 2.1 
31770000000000D-10,0.0000,6.35310000000000D- 10, 1.270620000000000- 09 
4,1. 76816150742336015, 1.23730971890897008, 1.45511293919343006,-3. 97 
57605 37 1 04600 DO 0, -1.50622 51 6081 086D- 02 ,4. 32 9 40 7 40 7326 48D- 04,-1. 85 35 
65607372189D-06,2.05286315303314D-09, 7. 36440814840596013, -5. 2584624 
73630271D06, 4.144787976812730-02/ 

P = PP 

T = TT 

C1=CC 

C2=C1*C1 

C3=C1*C2 

C4=C1*C3 

GO TO (1,7) , K 

1 GO TO (2, 3) , J 

2 N=*l 

GO TO 4 

3 All , 21-All, 1 1+AA/T 
A ( 2 , 2 ) = A( 2, l J-AB/T 

A ( 3,2 ) =A( 3, 1 ) +1.0D0/T 
A ( 5 , 2 ) =2. 0D0*A (3,2) 

N=2 

4 00 5 1=1,2 
M=2*N-2+I 

IF (M.EQ.2) GO TO 6 

F(l,M)=(A(l,N) + <A(2,N)+(A(3,N)+A44,N)*P)*P)*P}/[;i 
F(2,M)=(A(2,N)+(A(5,N)+A(6,N)*P)*P)/C2 
F(3,M)=(A(5,N)+A(7,N)*P)/C3 / 

F(4,H)=A(7,N)/C4 1 

5 P=U A 
RETURN 

6 F ( 1 ,2 ) =V1/C 1 
F ( 2 , 2) -V2/C2 
F (3,2) =V3/C3 
F ( 4 , 2 ) =F ( 4, l ) 

RETURN 

7 C5=C4*C1 
C6=C5*Cl 

GO TO (8,9) , J 

8 N = 1 

GO TO 10 

9 B(1,2)=B(1,1)+T *01 
B(2,2) =B( 2, 1)*T*D2 
8(3,2)=B(3,1)+T*D3 
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8(4,2>=8(4,1>+T*D4 
B(5,2)=BI5,1>*T*D5 
B<7,2) =B( 3,2)*2.0D0 
B(8,2)=B(4,2»*3.0D0 
B(9,2)=B(8,2>*2.000 
B< 10,21=8(5, 2)*4,ODO 
B(ll»2)=B(iO»2)*3,ODO 
B(12»2)=B(ll»2) *2.000 
N=2 

M=2*N-l 

G(l*M)=(B(i»N)MB(2»N) + (B(3»NI + (B(4,N) + (B(5»N)+B(6»N)*P)*PI*P)*P)* 

G(2»M)=(B(2,N)*(B(7,N)MB(B,N) + (B(10,NI+-B<13»N)*P>*P>*P>*P>/C2 
G(3,M)=(B(7,N)+(B(9,N)+(B( 11,N)+B(14,N)*P)*P)*P>/C3 

G ( 4 , M) = ( B ( 9 * N I ♦ < B(12»N)+B( 15, N ) *P I *P )/C4 
G ( 5 ,M>=(B(12,N)+B(16,N)*P)/C5 
G(6,M)=B< 16 » N >/C6 
M=M*l 

G(1»M)=B(1,N)/C1 

G(2,M) = B( 2 * N ) /C2 

G(3,MI=8(7,N»/C3 

G(4,MI=B(9,N)/C4 

G(5,HI=B(12,NI/C5 

G ( 6, M) =B( 16,N)/C6 

RETURN 

END 


C OECK METEMP 

DOUBLE PRECISION FUNCTION CP(T) 

DOUBLE PRECISION SI C 2 ) ,HI( 2 I * S 

DIMENSION M9.3.2) _ „ , 

DATA A/ 2. 5771104E-6, -2. 2 4078 IE- A* - A. 6776567E-4, 7 .7524692E-3 , 1 • 02 07 
1347E-3»-5.9827343E-2 * . 1053479, -6.7124682E-2 , 3. 01 59729, 3. 221388E- 7 , 
2-3.201ll57E-5,-7.7960945E-5,l.5504938E-3,2.5518368E-4,-1.9942448E- 
32,5.267395E-2,-6.7124682E-2,3.0159729,2.863456E-7,-2.8009762E-5,-6 

4.6823557E-5 ,1.2920782E-3,2.0414694E-4»-1.4956836E-2»3»5ll5967E-2»- 
53.3562341E-2,3.0159729,2.1771302E-6,-3.6924768E-5,1.2043213E-4,1.7 

6546467E-3,-l. 7244897E-2, 1. 88255 12E-2,. 4503988, -1.6311027, 4. 5834702 
7 t 2. 721 41 28E-7, -5. 2749669 E-6, 2.0QB0355E-5 , 3. 5092934E-4, -4. 31 12 242 E- 
83, 6, 275 1707 E-3,. 2251 994, -1,631 1027, 4, 5834702, 2.4190336E-7, -4. 61 559 
96E-6»1.7211733E-5,2, 9244112E-4, -3. 4489794E-3, 4.706378E-3 ,. 15013293 
$,-,81555135,4,5834702/ 

DATA SI,HI/ 18.66792402732497, 19.90897487890906, -2. 176323905587196, 
1-110.4372755187238/ 

K = 1 

1 N=1 

IFI T,GE,259«78828)N=2 
S=T/1, 3D2 
CP=A( i,K,NI 
DO 2 J=2, 8 

2 CP=CP*S+A< J,K,N) 

GO TO (3,4,51 , K 

3 CP=CP*S+A(9,1,N> 

RETURN 



entry csm 

K = 2 

GO TO l 

4 CP=CP*S+A(9,2,N)*DLOG<S)+Sl(N) 
RETURN 

ENTRY CH( T ) 

K = 3 

GO TO 1 

5 CP=T*(CP*S+A<9,3,N))+HI{N) 
RETURN 

END 


C DECK MELOG 


l 


LOGICAL FUNCTION LGFN{ P, T, J , Z ) 
COMMON /LDATA/ XKV, R , XMW, RC , D2»G 
DIMENSION Z<6) 

S^T/IOD.O 


J = 1 


LGFN=. TRUE. 

in^ i?^ T *^* lE7 * OR * P * LT,C>#l * OR * S * LT »- 69 * OR * s -GT.6>.1.0R.Zm.LE.O. 
10R.Z(2).LE.O.O.OR.Z( 3I.LE.0.0) RETURN 
IF (S.GT. 1.9077) GO TO 1 


PLOG^B. 30516+1 -2. 961— 0. 8/S ) / S 

IF (S.SE. 1.1883) PL0G=PLQG+0. 257* (S/1.1 88 3- 1.0)**!. 32 

IF ( P.GT.XKV*EXP(2.302585i*PLOG) ) RETURN 

J=0 


0 . 


LGFN=. FALSE. 

RETURN 

END 


C DECK MEDATA 

BLOCK DATA 
COMMON /LDATA/ R(61 

DATA R/1.0 T 518.2562* 16.04303* 8.745139E-4, 5.6, 1.33 33333/ 


C DECK METLG 

SUBROUTINE TLOGIC <P,T) 

DATA A1,A2, A3»A4, A5» A6» A7, A8, A9/53. 88758, 1. 82535 77, 0. 18723912 . 1 . 57 
10661E-5,-8.745l662E-4, 1 . 2470553E-4, 9.480861 7E-6, - 1 . 280319E -6, 4. 544 
26557E-B/ 

IF (T.GT.190.77) RETURN 
V=ALOG( P) 
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S = Al-M42*(A3-MA4 + (A5-MA6 + ( A7 + I A8+A9*V )*V ) *V I *V ) * V ) *V) *V ) *V 

IF (T.LT.S) T=S 

IF (T.LT.69.0) T = 69, 0 

RETURN 

END 


Natural Gas 

The routines for natural gas are described and presented in reference 3. The tem- 
perature is from 200 to 400 K, and the pressure range extends to 100x10^ pascals. The 
composition of natural gas is assumed to consist of alkanes containing from one to six 
carbon atoms and the dilutant gases Ng and CC^- The alkanes C^Hjq, C & Hj 2 , and 

can have more than one molecular configuration. For these cases, the assump- 
tion is made that each molecular configuration, for a given gas, is equally probable. 

The reference state for the enthalpy and entropy is that of the ideal gas at a temper- 

5 

ature of 200 K and a pressure of 1x10 pascals. 

Since natural gas is a mixture of many gases, an additional subroutine, whose deck 
name is NMCOMP, is necessary. This subroutine supplies composition-dependent con- 
stants to the other routines. For a given composition, this subroutine has to be refer- 
enced only once in the main program. This reference has to occur prior to a reference 
to either RWEDG or RNS. The reference to this subroutine is 

CALL BDATA(X) 

where X is an eight-element array whose elements are proportional to the mole frac- 
tions of the natural-gas components. The order in which these components appear in the 
array is as follows: CH^, CgHg, CgHg, C^Hjq, C & H 12 , C g H 14 , Ng, and COg. The 
natural-gas routine, as distinguished from the other routines, requires input of these 
values of X as additional data. 

Even though the card listing of these routines is in reference 3, this listing is 
repeated here for convenience. The prefix for the deck names is NM. The subroutine 
NMZETA is not independent; but the combination of NMZETA, NMPOLY, and NMCOMP 
is. The routine NMTEMP is not independent, but the combination of NMTEMP and 
NMCOMP is. 
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C DECK NMCOMP 


SUBROUTINE BDAT A (X) 

DIMENSION X<8), MOL18), XMOL ( 8 ) , SIS), HI8), CPI8,8), T(8), PI8), 
IRHOI8) 

COMMON /LOATA/ XKV, R , MW , RC , 02, G 

COMMON /PDATA/ FI 9 ) / ZOAT A/PC, TC, RHOC/TDATA/ A ( 8 , 3 > ,HI , S I 
REAL MOL, MW 

DATA MOL, XMOL/ 16. 043, 30, 07, *4. 097, 58. 12 A, 72. 151, 86. I?8, 28. 013, 44.0 
11 ,2 • 77527262, 3 •403528,3, 78 639175, 4. 06264748,4.27876115,4.45641492, 
23.33266869,3.78441688/ 

DATA CP/2 .7 998 255, .42 84998, -.2751805, 2. 5 82 171 IE- 2, 2.4165792 E-2 , -2. 
15163737E-3,-8.2465805E-4,l. 1523272E-4, -9. 8533835, 19. 657673, -10. 186 

2582. 1.8267443 .. 2463681, -.1202048, l. 0807487E-2, 0. 0,-16. 796807, 29. 08 
34569,-13.810883,2.2198327, .3655141, -.1532602, 1.296668E-2, 0.0, -1.81 
422946, 5. 6640979, -.907714,. 1435233,. 034644782, -.017195974, 1.7660626 
5E-3,0. 0,-3. 3598014, 7. 4196271, -.726671,. 045531 828, 4*0.0 ,-.5379224, 5 

6. 65 39353.. 3607859,-. 1684308, .015475231, 3*0.0, 2.501146,-9. 720581E-3 
7, 1. 036 056E-2,-4. 4372 58E-3, 6. 825596E-4, 3*0.0, 2. 5044684, -.5085567,. 4 
8840302,-3.73057lE-2,-2.522643E-2,6.1401476E-3,-4.1l66357E-4,0.0/ 

DATA S,H/-2. 42592233, -16. 722706, -24. 4685144, -7. 4352313, -9. 71 086973 
1,-9. 62405754, -1.20430845, -.54815092, -794. 255051, -224. 353146, 43. 254 
268,-792.772 573,-836. 398938,-1261.94398,-699.7098 35,-702.986595/ 
DATA P,T,RH0/4.626E6,4.894E6,4.257E6,3.722E6,3.299E6,3.149E6,3.398 
1E6,7.358E6, 190. 77, 30 5. 56, 369. 97, 416. 7, 454. 6, 499. 7, 126. 135, 304. 20,1 

262.5.203.2.220.5.224.4.235.0. 236.7.311.0.468.0, 

XX=0.0 

DO 1 N=l, 8 

1 XX=XX*XIN> 

DO 2 N=1 ,8 

2 F(N)=XIN)/XX 

FI9)=FIB)-F(7>/2.0+F12)+2.0*FI 3>+3.0*F( 4)*4.0*F(5 )+5.0*FI6) 

PC=0.0 
TC-0.0 
RHOC=0.0 
Sl=0.0 
H I =0.0 
MW=0. 0 
DO 3 N=1 » 8 

SI=$I+F(N)*(S(N) -XMOL I N ) ) 

HI = HI«-F(N)*HIN) 

MW=MW+FIN>*MOLIN) 

PC=PC+F(N)*P(N) 

TC=TC*F(N)*T IN) 

3 RHOC=RHOC+F ( N ) *RHOI N ) 

S I = S I +ALOGI MW) 

PC=P(i)/PC 
TC=T(i)/TC 
RHOC=RHO( 1 ) /RHOC 

DO 5 N=1 » 8 
NN=9-N 
A I NN, l ) =0.0 
DO 4 M=1 , 8 

4 A(NN,l)=A(NN,l )+F!M)*CP< N, M ) 

XN=N-i 

IF (XN.EQ.0.0) XN=1. 0 
AINN,2)=A(NN,1)/XN 

5 A(NN,3)=A(NN,1)/FL0ATIN) 
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R=S3l4. 41/MW 

RC=5.45105E-5*MW 

RETURN! 

END 


C DECK NMZETA 

SUBROUTINE ZETA (KK,PP,TT,Z) 

COMMON /VALUE/ F ( 4, 4 I , G t 6, 4 > 

COMMON /ZDAT A/ PC,TC,RHOC 
DIMENSION Z(6» 

DOUBLE PRECISION F,G , 61, B2,B3, B4, B5, Al, A2, A3 , A4, A5 ,E1 ,E2 ,PA , TA, THl 
1, TH2,TH3,TH4,D1,D2,03,D4,D5,F1,F2,UA,P, T, PI, P2 , U, T1 ,RC , FXPC , RR, E XP 
2B,ZB1,ZC1, AB1, AB2,AB3, AB4 r AB5, Z A, ZB, ZC, RBI , EXPB1 , S, SS , P$ I l , PSI 2 , PS 
313, PSI4,RC1» EXPCI, PSI5,PSI6,PSI7,PSIB 
DATA B1,B2,B3,B4,85, Al , A?, A3, A4, A5, El, E2, PA/ 4. 9147357499 1686D-C3,7 
1.376642234785500-06,-1. 145878430329230-07,5. 8951020951 11410-10,-5. 
2743822 81343532 D- 13, -2. 23 98 3199201862D0O, 1.343312537412700-03,2.759 
31 01 829 D6551D-05,-l. 655469770535420-07, 2. 34124562687064D-10, -4. 6002 
40000000000D-02, -2. 117700000000000-10,1.13318000000000002/ 

DATA TA,TH1,TH2»TH3, TH4, 01, 02, 03, 04, D5, F 1, F2 ,UA/ 1 . 4771 0550000 000 DO 
12,1.09934666473654D-14, 1.64873321284064007, 1.07243639762491008,3.6 
2644888 82455140-15, -3. 9776053710460DDOO, -1.506225160810860-02,4.329 
3407407 32648D-04, -1.853556073721890-06,2.052863153033140-09,-1.3787 
49330033000003,1. 3441 8460000000000,1 .4551 12 9391 9343D06/ 

K=KK 

P=PP*RHOC 

T=TT*TC 

P1=P+P A 

P2=P1*P1 

U=P2*P1 

T1=T+TA 

RC= I Fl*F2*P I /T 
EXPC=DEXPI RC J 
RB= ( E1+E2*U) 

EXPB=DEXP(RB*T1) 

ZB1=<TH1*P*P2*<U-TH2 >♦( TH3-U )*EXPB>/T 

ZCl = (Dl + ID2+( 03-M D4*D5*P )*P » *P ) *P > *P*EXPC/T 

IF (K.EQ.2I GO TO 1 

ABl=Bl*Al/T 

AB2=B2+A2/T 

AB3=B3+A3/T 

AB4=B4*A4/T 

AB5 =B5«-A5/T 

ZA=1.300+< ABl+t AB2*< AB3+( AB4+AB5+P ) *P ) *P )*P J*P 
Z ( 1 ) =Z A+-ZB1 + ZC1 
IF (K.EQ.O) RETURN 

ZA=1.DOO+(2.0DO*AB1M 3. 0D0*AB2M 4. QD0*AB3+< 5. QD0*AB4+6. 0D0*AB5*P > * 
1P)*P)*P)*P 

ZB=ZBl*(2.0DO*(1.0DO+P/Pl>+3.0DO*P*P2*(E2*Tl+1.0DO/<U-TH3H-l.OOO/< 
1U-TH2) I J 

ZC= I F2*P*ZCl+£XPC*<2.OD0*Dl + (3.OD0*D2M4.ODO*O3+< 5. 000*04+6. 000*D5 
1*P)*P)*P)*P)*P|/T 
Z ( 3 ) =Z A+Z8+ZC 
IF (K.EQ.l) RETURN 
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1 


RB1=E1«-E2*UA 
EXPt5l = DEXP( RR1 *T1 ) 

ZA=l.ODO+(Bl+(B2+{83+(B4+B5*P)*P)*P>*P)*P 

ZB=RB*T*ZB1 

ZC=-RC*ZC1 

Z ( 2 ) =Z A+ZB+ZC 

S=E2*T 1 

SS=F2/T 

CALL POLY ( 1, 1,U,T,S) 

CALL POLY (1,2,P,T,SS) 

ZA=< BH-(82/2.0D0+( B3/3.0D0M B4/4,ODO+B5*P/5.QOO) *P)*P)*P)*P 
P$ll=F(l,l)-FC2,l)*F<3,l)-F(4,l) 

P$I2=F<1,2)-F<2,2)+F(3,2)-F(4,2) 

ZB=TH4*(P$I 1*EXPB-PS I2*EXPB1 ) 

PSI3=G<1,1)-G(2,1)*SC3,1)-G<4,1)+G(5,1)-G<6,1) 

PSI4=G(1,2)-G( 2,2)+G( 3,2 )-G< 4,2)+G< 5,2)-G< 6,2) 

RCi-Fl/T 
EXPC1=DEXP( RC1 ) 

ZC=< PSI4*EXPC1-PSI3*EXPC)/T**2 
Z(4) =ZA+ZB+ZC 
IF (K.EQ.2) RETURN 
CALL POLY <2,i,U,T,S) 

CALL POLY (2,2,P,T,SS> 

PSI5=F(1, 3>-F(2,3)+F(3,3)-F(4,3) 

P$I6=F(i,4)-F(2 f 4>+F(3,4)-F(4,4) 
PSI7=G<l,3)-G<2,3)+G(3,3)-G<4 # 3l+G(5, 3>-G<6,3> 
PSI0=G(1,4)-G(2,4)+G( 3,4)-G< 4,4)+G< 5,4)-G( 6,4) 

ZA = -( Al + < A2/2.000+< A3/3.0DQ+( A4/4. ODO+A5*P/5,ODO ) *P > *P ) *P ) *P/T 
ZB=TH4*(PSI5*EXPB-PS I 6*EXPB1 ) 

ZC=( PSIB*EXPC1-PSI7*EXPC)/T**2 
Z ( 5 ) =Z A+ZB+ZC 

ZB=TH4*T*(EXPB*(RB*PSI1-(F< 1 t 1)-2.0D0*F( 2,1 ) +3. 3DO*F < 3 , 1 )-4,ODO*F( 
14,1 ) )/Tl)-EXPBl*(RBl*PSI2-(F(l,2)-2.0D0*F( 2, 2) +3. ODO*F ( 3 , 2 ) -4. OOO* 
2F ( 4 , 2 ) )/Tl) ) 

ZC=< EXPC*( < 2.0DQ«-RC)*PSI3-G< 1, 1>+2.0DD*G< 2, 1 )-3. ODO*G< 3, 1 )+4. ODO*G 
1(4,1 )-5.*G< 5,1)+6.0DO*G( 6, 1 ) l-EXPC 1* ( ( 2 . ODO+RC 1 ) *PS I 4-G ( 1,2) 4-2. ODD 
2*G<2,2>-3.0DO*G(3,2>+4.00D*G(4,2)-5.0DO*G(5,2)+6.0DO*G(6,2) ) )/T**2 
Z (6 ) =ZB+ZC 
RETURN 
END 


C DECK NMPOLY 


SUBROUTINE POLY ( J, K , PP, TT , CC ) 

COMMON /VALUE/ F ( 4, 4 ) , G( 6, 4 ) 

DIMENSION A ( 7, 2 ) , B(16,2) 

DOUBLE PRECISION PP, TT , CC, F, G, A , B, AA, AB, UA , D1 , D2 ,D3»D4,D5,P,T,C1,C 
12,C3,C4,C5,C6,V1,V2,V3 

DATA B/5, 48429563564224D03, 1 .54210957 155370001 ,- 6 . 17182239272157D- 
101,3.13762297202746D-03,-5.32199549075239D-06,2. 75942703621 458D-09 
2,-1 . 2343644 78 54431 DDO, 9, 412868 916082 3 8D-03, 1.8825737832164BD-02,-2 
a q 707 J > ? J q * 386 3945 88902 86D— 05, -1.27727891778057D— 04, 1, 

5li^cTc2I 2 -2"2 8 * 5 * 5 l 805 4O72429 17 D-O8, 1.655656221728750-07,3.311 
61 + 0DC0 * 2 * 7594270 3621 4580-09 , 6*0 , DDOO , 1 . 37971 35 18 

61 0729D -OB, 5. 51B854072 429 17D-08, 1.655656221728750-07,3,311312443457 
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7500-07/ 

DATA 4,AA, A8, J A, D1 , D2 , 03, 04, D5 , V1,V2, V3/8. 133896 56644895 D1 3, -5. 31 7 

142860549802006, 1.979 949208 266470-02,2.117700000000000-1 0t 3. 95 9898 A 

21 65 32 93D-02 ,6.353100000000000-10, 1.270620000000000-09,3*0. 0D00,2 . 1 
31 770 OODOOOODOO-IO, 0.0000,6.353100000000000- 10, 1.270620000030000-09 
4, 1. 768 161 50742336D15, 1.2 3730971890897D08, 1.455 11 293919343006, -3. 97 
57605371046DODOO,-1.50622516081086D-02,4. 32940740732648D-04,-1. 8535 
6560737 2189D-06, 2. 052 8631 53033140-09, 7. 3644081484059601 3, -5. 2584624 
73630271006,4.144787976812730-02/ 

P = PP 

T = TT 

C1 = CC 

C2=C1*C1 

C3=C1*C2 

C4=C1*C3 

GO TO ( 1 » 7 ) , K 

1 GO TO (2,3) , J 

2 N-l 

GO TO 4 

3 A(1*2)=A(1»1)+AA/T 
A(2,2)=A(2,1)-AB/T 
A(3,2)=A(3»1)«-1.0D0/T 
A(5,2)=2.0D0*A(3,2) 

N = 2 

4 DO 5 1=1,2 
M=2*N-2+I 

IF (M.EQ.2) GO TO 6 

Ftl*M)«(A(l'N>MA(2'N) + (A< 3, N)*A< 4, N ) *P ) *P ) *P > /Cl 
F(2,M)=<A(2*N)+(A(5,N)+A(6»N)*P)*P)/C2 
F(3,M)=(A(5,N)+A(7,N)*P)/C3 
F ( 4 , M) =A(7,N)/C4 

5 P=UA 
RETURN 

6 F < 1 ,2) =V1/C1 
F (2,2) =V2/C2 
F ( 3 , 2 ) =V3/C3 
F ( 4 , 2 ) =F ( 4, 1 ) 

RETURN 

7 C5=C4*C1 
C6=C5*C1 

GO TO (8,9) , J 

8 N= 1 

GO TO 10 

9 B(1,2)=B1 1, 1 )+T*Ol 
B ( 2,2) =B( 2,1 >+T*02 
B ( 3 , 2 ) =8( 3, 1)+T*D3 
B(4,2)=B(4,1)+T*D4 
B(5,2)=B(5,1)*T*D5 
B(7,2)=B(3,2)*2.0DO 
B(8,2)=BC4,2)*3.000 
B(9,2)=B(8,2)*2.0D0 
B( 10,2)=B(5,2)*4.0DO 
B(il,2)=B(10,2)*3.000 
B(12,2)-B(ll,2)*2.000 
N = 2 

10 M=2*N-l 

GU,«) = (B(l,NM-(B(2,N) + ( B( 3, N )♦( B( 4, N ) + ( B( 5, N ) +6 ( 6, N) *P ) *P) *P ) *P ) * 
IP) /Cl 

G(2,MI=(B(2,N)+(B(7,N)+(B(8,N)+(B( 10,N)+B( 13,N)*P)*P)*P)*P)/C2 
G(3,M)=(B(7,N)+(B(9,N)+(B(11,N)+B(14,N)*P)*P)*P)/C3 
G(4,M)=(B(9,N)+(B(12, N)+B( 15, N )*P )*P )/C4 
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G(5,M)={B<12,N)+B(16,N)*P)/C5 
G(6,M) = B ( 1 6 , N ) / C6 
M = M+1 

G(l t M)*BCl»N)/Cl 

G(2,M)=B(2,N)/C2 

G(3»M)=8(7»N)/C3 

G(4,M)=B(9,N)/C4 

G ( 5 » M) =B( 12,N)/C5 

G ( 6 » M) = 8 ( 16»N)/C6 

RETURN 

END 


C DECK NMTEMP 

DOUBLE PRECISION FUNCTION CP(T) 
DOUBLE PRECISION S 
COMMON /TO AT A/A(B,3),HI,SI 
K = 1 

1 S=T/1.DD2 
CP=A(1,K) 

DO 2 N=2 * 7 

2 CP=CP*S*A ( N»K ) 

GO TO (3,4,5), K 

3 CP = CP*S«-A(8,1 ) 

RETURN 

ENTRY CS(T> 

K=2 

GO TO 1 

4 CP=CP*S + A(8,2)*DLOG(S)«-Sl 
RETURN 

ENTRY CH(T) 

K = 3 

GO TO I 

5 CP=T*(CP*S+A(8,3) )+HI 
RETURN 

ENO 


C DECK NMLOG 

LOGICAL FUNCTION LGFN ( P, T, J , l ) 

COMMON /ZDATA/ PC»TC,RHOC 
COMMON /LOATA/ XKV, R, XMW,RC, 02, G 
DIMENSION Z ( 6 ) 

S=T/ 103.0 
J=1 

LGFN=. TRUE. 

IF (P,3T.1.1E7.0R,P,LT.O,1.0R.S,LT»1.9.0R,S,GT.4.1,OR.Z(1 
IOR.Z(2).LE.O.O.OR.Z(3).LE.O.O) RETURN 
S=S*TC 

IF ( S.GT. 1.9077) GO TO 1 

PLOG=8.30516+(-2. 961 -0.8/S)/S+.257*(S/l, 1883-1.0)**!. 32 


) . LE. 0.0. 
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IF (P*PC.GT.XKV*EXP(2.3025851*PL0G)) RETURN 
1 J=0 

LGFN=. FALSE. 

RETURN 

ENO 


C DECK NMDATA 

BLOCK DATA 
COMMON /LDATA/ RI6) 

DATA R<1» ,RI5) ,RI6)/ 1.0, 5.6, 1.3333333/ 
END 


C DECK NMTLG 

SUBROUTINE TLOGIC IP,T) 

COMMON /ZDATA/ PC*TC,RHOC 
DIMENSION A ( 9 ) 

DATA A/4.5446557E-8»-1.280319E-6,9.4808617E-6, 1. 24705 5 3E -4, -8. 7451 
1662E-4.1.570661E-5,. 1872 3912, 1.8253577, 53. 88758/ 

PP=P*PC 

TT=T*TC 

IF (TT.GT. 190.77) RETURN 
V=ALOG<PP) 

S = 0.0 

DO 1 N=1 , 9 
l S=S*V* AIN) 

IF (TT.LT.S) T=S/TC 

IF(T.LT.190.0)T=190.0 

RETURN 

END 


81 



APPENDIX E 


SAMPLE CALCULATIONS 

It is desired to determine the properties of a supersonic air stream. These proper- 
ties are toe velocity; the Mach number; and the pressure, density, and temperature 
under both total and static conditions. Two cases are considered: one uses the routines 
in appendix B, and the other uses the routines in appendix C. For both cases, the input 
and output variables are in U.S. customary units rather than SI units. That is, pressure 
is in psia, temperature in degrees F, density in lb/ft , and velocity in ft/sec. Appropri- 
ate conversion factors are included in the main computer programs. 


CASE I 


For case I, the measurements are 

(1) The pressure indicated by a total-head tube. This is the total pressure down- 
stream of the normal shock. 

(2) The stagnation temperature, as indicated by thermocouple probe whose recovery 
factor is unity. This is the total temperature downstream of the normal shock. 

(3) The pressure indicated by a wedge -type static -pressure probe whose total wedge 
angle is 15°. This is the pressure behind the oblique shock attached to the wedge. 

These measurements are such that the routines in appendix B and the air routines 
in appendix D apply. However, in order to activate these routines, a main computer 
program must be written. A typical main program could be 

$IBFTC MAINW 

COMMON/OUTW/TP(8,5), VM(10), CONV(12), Z(6, 5),KODE(ll) 

REAL Ml 

NA ME LIST/OUT/ P2T, T2T, P3,DELTA2, Ul, Ml, PIT, RHOIT, TIT, PI, RHOl, 

1T1, KODE 
1 RE AD (5, OUT) 

P=P2T*6894. 73 
PP=P3*6894. 73 
T=(T2T+459. 67/1. 8 
DELTA=DELTA2/2. 0 
CALL RWEDG(P, T, PP, DELTA) 

Ul=VM(3)*3. 28083 
Ml=VM(l) 

PlT=TP(l, l)/6894. 73 
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Pl=TP(l,2)/6894. 73 
RH01T=TP(2, 1)*0. 0624284 
RHOl=TP(2, 2)*0. 0624284 
TlT=TP(3, 1)*1. 8-459.67 
Tl=TP(3, 2)*1. 8-459. 67 
WRITE (6, OUT) 

GO TO 1 
END 


The decks that are loaded are as follows; MAINW, RGWED, RGDFW, PGWED, 
WDATA, AIZETA, AITEMP, AILOG, AIDATA, and the necessary matrix inversion rou- 
tine deck called MINV which is supplied by IBM in their Scientific Subroutine Package. 

The input for a typical case is 

P2T Pressure indicated by the total head tube, 1027 psia 

T2T Temperature indicated by a stagnation thermocouple, 83. 4° F 

P3 Pressure indicated by a wedge -type static -pressure probe, 354 psia 

DELTA 2 Wedge angle, 15° 

The output for this case is 
Ul Velocity, 1517 ft/sec 

Ml Mach number, 1. 702 

PIT Total pressure, 1186 psia 

RHOIT Total density, 5.902 1b/ft^ 

TIT Total temperature, 86. 7° F 

Pi Pressure, 241. 1 psia 

RHOl Density, 2.005 lb/ft^ 

T1 Temperature, -119. 1° F 


CASE n 

For case n, the flow issues from a supersonic nozzle attached to a plenum . The 
assumption is made that the flow is isentropic and one dimensional from the plenum to 
the nozzle exit. The measurements are 

(1) The pressure in the plenum. This is the total pressure. 

(2) The temperature in the plenum. This is the total temperature. 
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(3) The pressure indicated by a total-head tube. This is the total pressure down- 
stream of the normal shock. 

The measurements are such that the routines in appendix C and the air routines in 
appendix D apply. However, in order to activate these routines, a main computer pro- 
gram must be written. A typical main program, which is similar to that in case I, 
could be 

$IBFTC MAINNS 

COMMON/OUTNS/TP (8 , 4) , VM(6) , CONV(9) , Z (6 , 4) , KODE (8) 

REAL Ml 

NAMELIST/OUT/PIT, TIT, P2T, Ul, Ml, RHOIT, Pi, RHOl, Tl,KODE 
1 READ (5, OUT) 

P=P1T*6894. 73 
PP=P2T*6894. 73 
T=(TlT+459. 67)/l. 8 
CALL RNS(P, T, PP) 

U1=VM(3)*3. 28083 
Ml=VM(l) 

P1=TP(1, 2)/6894. 73 
RH01T=TP(2, 1)*0. 0624284 
RHOl=TP(2, 2)*0. 0624284 
T1=TP(3,2)*1. 8-459. 67 
WRITE (6, OUT) 

GO TO 1 
END 

The decks that are loaded are as follows: MAINNS, RGNS, RGDFNS, PGNS, 
NSDATA, AIZETA, AITEMP, AILOG, AIDATA, and the necessary matrix inversion 
deck called MINV supplied by IBM in their Scientific Subroutine Package. 

The input for this case is 

PIT Plenum pressure, 1186 psia 

TIT Plenum temperature 86. 7° F 

P2T Pressure indicated by a total head tube, 1027 psia 

The values of PIT and TIT are given the same as those calculated in the previous 
case, also the value of P2T is the same as in the previous case. 

The output for this case is 

Ul Velocity, 1517 ft/sec 
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Ml 

Mach number, 1. 702 

RHOIT 

Total density, 5.902 lb/ft 1 

PI 

Pressure, 241.1psia 

RHOl 

Density, 2. 005 lb/ft 3 

Tl 

Temperature, -119. 1° F 


The results for this case are the same as in the previous case, as would be expected 
from the input. This similarity constitutes a good check on the validity of the routines 
in appendixes B and C . 
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